Hello,

This is my first attempt to implement a modeling simulation of neurons using Brian2. I would like to implement the simulation of the effect of ultrasound stimulation.

However, when the code is executed, the following error exists: the variable ‘U’, cannot be separated into linear components.

```
from brian2 import *
from scipy import integrate
from sympy import Integral
import numpy as np
newton = kgram * meter / second ** 2
Pa = newton / meter ** 2
T = 309.15 * kelvin
cm0 = 1 * uF / cm2
a = 32 * nmeter
delta = 1.26
A = 1 * Pa
f = 1000 / ms
zeR = 0 * nmeter
oneR = 1 * nmeter
z0 = 0 * nmeter
u0 = 0 * nmeter/ms
Rg = 8.314 * meter ** 3 / mole / kelvin
delta0 = 2.0 * nmeter
Delta = 1.26 * nmeter
Delta_ = 1.4 * nmeter
pDelta = 1.0e5 * Pa
m = 5.0
n = 3.3
rhoL = 1075 * kgram / meter ** 3
muL = 7e-4 * Pa * second
muS = 0.035 * Pa * second
kA = 0.24 * newton / meter
alpha = 7.56 * Pa * second
C0 = 0.62 * mole / meter ** 3
kH = 1.613e5 * Pa * meter ** 3 / mole
P0 = 1e5 * Pa
Dgl = 3.68e-9 * meter ** 2 / second
xi = 0.5 * nmeter
epsilon0 = 8.854e-12 * farad / meter
epsilonR = 1.0
rel_Zmin = -0.49
Z_qs = 0.001 * nmeter
V_T = -56.2 * mV
gNa = 56 * mS / cm2
vNa = 50 * mV
gK = 6 * mS / cm2
vK = -90 * mV
gM = 0.075 * mS / cm2
vM = -90 * mV
gLeak = 0.075 * mS / cm2
vLeak = -70.3 * mV
def localDeflection(r, Z, R):
if np.abs(Z) == 0.0:
return 0.0
else:
return np.sign(Z) * (np.sqrt(R ** 2 - r ** 2) - R + Z)
def PMlocal(r, Z, R):
z = localDeflection(r, Z, R)
relgap = (2 * z + Delta) / Delta_
return pDelta * ((1 / relgap) ** m - (1 / relgap) ** n)
def PMavg(Z, R, S):
fTotal, _ = integrate.quad(lambda r, Z, R: 2 * np.pi * r * PMlocal(r, Z, R),0, a, args=(Z, R))
# for num in range(1000):
# r1 = a / 1000 * num
# r2 = a / 1000 * (num + 1)
# j1 = 2 * np.pi * r1 * PMlocal(r1, Z, R)
# j2 = 2 * np.pi * r2 * PMlocal(r2, Z, R)
# if not num:
# fTotal = (j1 + j2) * a/1000/2
# else:
# fTotal += (j1 + j2) * a/1000/2
return fTotal / S
def surface(Z):
return np.pi * (a**2+Z**2)
def gasFlux(Z, P):
''' Gas molar flux through the sonophore boundary layers
:param Z: leaflet apex deflection (m)
:param P: internal gas pressure (Pa)
:return: gas molar flux (mol/s)
'''
dC = C0 - P / kH
return 2 * surface(Z) * Dgl * dC / xi
PMavg = Function(PMavg, arg_units=[meter, meter, meter**2], return_unit=Pa)
gasFlux = Function(gasFlux, arg_units=[meter, Pa], return_unit=mole/second)
eps = """
dZ/dt = U : meter
dU/dt = -3/2/R*U**2 + 1/rhoL/abs(R)*Ptot : meter/second
Ptot = Pin+Pm+Pec-P0+PA-Ps-4/abs(R)*U*(3*delta0*muS/abs(R)+muL) : kgram / meter / second ** 2
Pin = int(Z < Z_qs)*(P0*pi*a**2*Delta/Va) + int(Z >= Z_qs)*(P0*ng*Rg*T/Va) : kgram / meter / second ** 2
dng/dt = gasFlux(Z, Pin): mole
Pm = PMavg(Z, R, S) : kgram / meter / second ** 2
Pec = -pi*a**2/S * (c_m*v)**2/(2*epsilon0*epsilonR) : kgram / meter / second ** 2
Ps = 2*kA*Z**3/(a**2*(a**2+Z**2)) : kgram / meter / second ** 2
c_m = cm0*Delta/a**2*(Z+(a**2-Z**2-Z*Delta)/(2*Z)*log((2*Z+Delta)/Delta)) : farad / metre ** 2
dv/dt = -(gNa * m ** 3 * h*(v-vNa)+gK * n ** 4 *(v-vK)+gM * p*(v-vK)+gLeak*(v-vLeak))/c_m : volt
dm/dt = alpham * (1-m) - betam * m : 1
dn/dt = alphan * (1-n) - betan * n : 1
dh/dt = alphah * (1-h) - betah * h : 1
dp/dt = alphap * (1-p) - betap * m : 1
alpham = -(0.32/mV)*(v-V_T-13*mV)/(exp(-(v-V_T-13*mV)/(4*mV))-1)/ms : Hz
betam = (0.28/mV)*(v-V_T-40*mV)/(exp((v-V_T-40*mV)/(5*mV))-1)/ms : Hz
alphah = 0.128 * exp(-(v-V_T-17*mV)/(18*mV))/ms : Hz
betah = 4/(exp(-(-v-V_T-40*mV) / (5*mV)) + 1)/ms : Hz
alphan = -(0.032/mV) * (v-V_T-15*mV)/(exp((-v-V_T-15*mV)/(5*mV))-1)/ms : Hz
betan = 0.5*exp(-(v-V_T-10*mV)/(40*mV))/ms : Hz
alphap = 1/(1+exp(-(v+35*mV)/(10*mV)))/ms : Hz
betap = 608/(3.3*exp((v+35*mV)/(20*mV))+exp(-(v+35*mV)/(20*mV)))/ms : Hz
R = (a ** 2 + Z ** 2)/(2 * Z) : meter
S = pi * (a**2 + Z**2) : meter ** 2
Va = pi * a**2 * Delta * (1+Z/(3*Delta)*(Z**2/a**2+3)) : meter ** 3
r : meter
PA : kgram / meter / second ** 2
"""
legs = NeuronGroup(1, model=eps, method='exponential_euler')
trace = StateMonitor(legs, 'c_m', record=True)
run(100*ms)
plot(trace.t/ms, trace.c_m[0])
```