Simulation of ultrasound stimulation using Brian2 fails to compile nonlinear formulas

Support
1

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])
2

Hi @lay. Interesting project ! Regarding the error, the full message states:

…
ValueError: The expression ‘…’, defining the variable ‘U’, could not be separated into linear components.
…
brian2.stateupdaters.base.UnsupportedEquationsException: Can only solve conditionally linear systems with this state updater.
…

This means that for these equations, you cannot use the exponential_euler algorithm, since your equation for dU/dt depends on U^2. Instead, you could choose e.g. method='rk2' or method='rk4' to use a Runge-Kutta method.

This will unfortunately not make things run right away, though. First, since you are defining functions in Python, you cannot use the Cython/C++ code generation and should state

prefs.codegen.target = 'numpy'

Then, you need to handle unit issues in the functions. The numpy functions will work fine with the units (e.g. np.sqrt(R ** 2 - r ** 2) would give you a value in meter), but the integrate.quad function will not. It will not touch the arguments Z and R (so they will arrive as quantities with units in the function PMlocal), but r will not have any units attached. I think the easiest way is to use

fTotal, _ = integrate.quad(lambda r, Z, R: 2 * np.pi * r * PMlocal(r, Z, R), 0, a/nmeter, args=(Z, R))

i.e. give unitless values for r to the function, and then reattach units in the function PMlocal:

z = localDeflection(r*nmeter, Z, R)

The value of fTotal will also be without units, so you need to use something like

return fTotal*newton / S

Note that you do not actually define a gasFlux function, you could directly include gasFlux and surface in the equations, e.g. as:

surface = pi*(a**2+Z**2) : meter ** 2
dC = C0-Pin/kH : mole / meter ** 3
dng/dt = 2 * surface * Dgl * dC / x_i: mole

Note that I replaced xi by x_i, since the xi variable has a special meaning in Brian (noise term).

With all these changes, your simulation should run, but it will not give any meaningful results yet. The reason is that the initial value for Z is zero, which in turn means that R divides by zero, and the value of a lot of variables turns into inf or NaN. I don’t know these equations well enough to recommend any specific course of action, maybe setting Z to an initial non-zero value is enough?

Hope that gives you some pointers!

3

The part I used is Equation (1)-Equation (11) in the figure. In the subsequent run it does appear that the result is Nan, may I know in what way I need to initialize the Z value.

bls
bls626×833 198 KB

4

Hi @lay. I am afraid I am unfamiliar with this type of equation and do not know how they are commonly used. But from the equation, you see that R(Z) is a function like this

Figure_1
Figure_1640×480 9.06 KB

For Z equal to 0, i.e. its initial value, it is undefined. You’ll need to either start it at a non-zero value or change the equations so that they are defined for all values of Z.