# Implementing second-order differential equations in Brian2

**URL:** https://brian.discourse.group/t/implementing-second-order-differential-equations-in-brian2/450
**Category:** Science, Projects, and Showcases
**Tags:** equations
**Created:** [11 August 2021 07:48 UTC](https://brian.discourse.group/t/implementing-second-order-differential-equations-in-brian2/450 "2021-08-11T07:48:57Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Kim](https://avatars.discourse-cdn.com/v4/letter/k/9f8e36/32.png) [@Kim](https://brian.discourse.group/u/Kim)
#### Post date: [11 August 2021 07:48 UTC](https://brian.discourse.group/t/implementing-second-order-differential-equations-in-brian2/450/1 "2021-08-11T07:48:57Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![adam](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/adam/32/197_2.png) [@adam](https://brian.discourse.group/u/adam)
#### Post date: [12 August 2021 15:04 UTC](https://brian.discourse.group/t/implementing-second-order-differential-equations-in-brian2/450/2 "2021-08-12T15:04:35Z")

</div>

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](https://brian2.readthedocs.io/en/stable/user/converting_from_integrated_form.html)

from [Equations](https://brian2.readthedocs.io/en/stable/user/equations.html)

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](https://brian2.readthedocs.io/en/stable/advanced/namespaces.html)

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

```