Implementing second-order differential equations in Brian2

Science, Projects, and Showcases
1

Hello, I have some question about implementing the second-order differential equations to the Brian2.

The question I want to ask you while using the Brian2 is:
Is there any way that you can implement the second-order differential equation for the model of the synapse in Brian2?
Or, at least get the t in R(t), without getting inverse function. The inverse function of R(t) could not be calculated from the MATLAB or Wolframalpha or evenwith any other methods I tried. Thus, I found the pynverse package but the Brian2 could not resolve it.

This is the equation I’ve wanted to use for the synapse model, but I figured Brian2 can only handle the first-order differential equations. The second line is the original form of the equation, the solution of the differential equation.

𝑅^β€²β€²+16.1424𝑅^β€²+16.2660𝑅=16.2660

𝑅=βˆ’0.1126 𝑒^((βˆ’π‘₯/0.926) )βˆ’0.282 𝑒^((βˆ’π‘₯/0.06639) )+1

So, I used the internal variable z to make this equation into two first-order differential equations, that look like these:

𝑑𝑅 / 𝑑𝑑 = 𝑧 βˆ’ 16.1424 𝑅
𝑑𝑧 / 𝑑𝑑 = 16.2660 βˆ’ 16.2660 𝑅

where,

𝑍 = 𝑅′ + 16.1424 𝑅
= βˆ’ 1.6960375 𝑒^(βˆ’π‘₯ / 0.926) βˆ’ 0.3045318𝑒^(βˆ’π‘₯ / 0.06639) + 16.1424

Until here, it was working fine, no error was happening. However, the signal R changes abruptly when the presynaptic spike happens. It decreases to 0.6R, which makes signal z also have to change abruptly. When I was trying to get the value or equation for updating z, I found the problem.

Since z is the internal variable to calculate R, z and R should be calculated using the same t. But as you know, if I put these two equations like that, and also shown below, Brian2 will consider z and R as a separate signal and will use different t on calculating them.

Plus, to get the z value when the R becomes 0.6R, I need to calculate R’ where t= (0.6R)^-1. I was trying this way and found pynverse package that calculates the inverse function’s value on the specific point. The problem here is that when I use pynverse package and trying to put that value like β€˜z = invR(R)’, Brian2 rises an error saying that invR could not be resolved.

The part of the code that I am working on is like this. And this rises the error.

#######################################################
from pynverse import inversefunc

funcR = lambda x : - 0.1126 * np.exp(- x / 0.926)  - 0.282 * np.exp(- x / 0.06639) + 1
invR = inversefunc(funcR, domain=-1, open_domain=True)

U = 0.6054
taur1 =  1000 * ms
taur2 =  1000 * ms

SpGe_L = Synapses(in_L, L, 
             model = '''
                      dR/dt =            ( z - 16.1424 * R ) / taur1 : 1 (clock-driven)
                      dz/dt  = ( 16.2660 - 16.2660 * R ) / taur2 : 1 (clock-driven)
                     ''', 
             on_pre = ''' 
                      x_e_post += w1 * e * R
                      R *= U
                      z = invR(R)
                      ''',
              name = 'in_L', method = 'euler') 

SpGe_L.connect(j = 'i') 
SpGe_L.z  = 16.1424 
SpGe_L.R  = 1

#######################################################

(I guess I cannot upload the figures here)

Whenever R changes abruptly, z should be recalculated with the same t that is used for new R(t). I was trying to find this t with invR(R), but failed.

In this situation, I could think of two solutions:

  1. Implement the second-order differential equation with R’’ directly into the Brian2, if possible (this is the best way)

  2. Find a way to use pynverse package in the model equation of the Brian2, like I did β€˜z = invR(R)’.

Any advice will be helpful for me.

2

Hi JKim,
(caveat, i’m neither a brian expert, nor a numerical integration expert)

From looking at the docs, I don’t think there’s a way to implement the model as a second-order differential equation apart from breaking it into first order differential equations as you’ve done.

Brian models are defined by systems of first order ordinary differential equations, but you might see the integrated form of synapses in some textbooks and papers. See Converting from integrated form to ODEs

from Equations

If using invR( ) ends up being the right way to go, I think this might be a simple scope issue. I think one way to solve this is to specify invR in the namespace of the simulation when you run it

docs for Namespaces

net = Network()
net.run(sim_duration, namespace={'invR':invR})