# Simple model using double exponentials


```python
from sympy import *
```


```python
from IPython.display import display, Markdown
```


```python
init_printing()
```


```python
t, P, e_r, e_d, delta_e, rho_e, g_e, i_r, i_d, delta_i, rho_i, g_i, b = symbols('t P \\tau_{er} \\tau_{ed} \\delta_e \\rho_e \\bar{g}_e \\tau_{ir} \\tau_{id} \\delta_i \\rho_i \\bar{g}_i \\beta', positive=True, real=True)
```


```python
SymbolDict = {t: "Time (ms)", P: "Proportion of $g_i/g_e$", e_r: "Excitatory Rise (ms)", e_d: "Excitatory Fall (ms)", delta_e: "Excitatory onset time (ms)", rho_e: "Excitatory $tau$ ratio (fall/rise)", g_e: "Excitatory max conductance", i_r : "Inhibitory Rise (ms)", i_d : "Inhibitory Fall(ms)", delta_i:"Inhibitory onset time(ms)", rho_i: "Inhibitory $tau$ ratio (fall/rise)", g_i: "Inhibitory max conductance", b: "Inhibitory/Excitatory $tau$ rise ratio"}
```


```python
estimateDict = { P: (1.9,2.1), e_r: (1.5,5), e_d: (8,20), delta_e: (0,0), rho_e: (2,7), g_e: (0.02,0.25), i_r : (1.5,5), i_d : (14,60), delta_i:(3,8), rho_i: (5,20), g_i:(0.04,0.5) , b:(0.5,5) }
```


```python
averageEstimateDict = {key: pow(value[0]*value[1],0.5) for key,value in estimateDict.items()}
```


```python
#averageEstimateDict = {key: ((value[0]+value[1])/2.) for key,value in estimateDict.items()}
```


```python
print "| Variable  |  Meaning |  Range |"
print "|---|---|---|"
print "|$t$|Time (ms)|0-100|"
for i in [P, e_r, e_d, delta_e, rho_e, g_e, i_r, i_d, delta_i, rho_i, g_i, b]:
    print "|${}$|{}|{}-{}|".format(i, SymbolDict[i], estimateDict[i][0], estimateDict[i][1])
```

    | Variable  |  Meaning |  Range |
    |---|---|---|
    |$t$|Time (ms)|0-100|
    |$P$|Proportion of $g_i/g_e$|1.9-2.1|
    |$\tau_{er}$|Excitatory Rise (ms)|1.5-5|
    |$\tau_{ed}$|Excitatory Fall (ms)|8-20|
    |$\delta_e$|Excitatory onset time (ms)|0-0|
    |$\rho_e$|Excitatory $tau$ ratio (fall/rise)|2-7|
    |$\bar{g}_e$|Excitatory max conductance|0.02-0.25|
    |$\tau_{ir}$|Inhibitory Rise (ms)|1.5-5|
    |$\tau_{id}$|Inhibitory Fall(ms)|14-60|
    |$\delta_i$|Inhibitory onset time(ms)|3-8|
    |$\rho_i$|Inhibitory $tau$ ratio (fall/rise)|5-20|
    |$\bar{g}_i$|Inhibitory max conductance|0.04-0.5|
    |$\beta$|Inhibitory/Excitatory $tau$ rise ratio|0.5-5|


| Variable  |  Meaning |  Range |
|---|---|---|
|$t$|Time (ms)|0-100|
|$P$|Proportion of $g_i/g_e$|1.9-2.1|
|$\tau_{er}$|Excitatory Rise (ms)|1.5-5|
|$\tau_{ed}$|Excitatory Fall (ms)|8-20|
|$\delta_e$|Excitatory onset time (ms)|0-0|
|$\rho_e$|Excitatory $tau$ ratio (fall/rise)|2-7|
|$\bar{g}_e$|Excitatory max conductance|0.02-0.25|
|$\tau_{ir}$|Inhibitory Rise (ms)|1.5-5|
|$\tau_{id}$|Inhibitory Fall(ms)|14-60|
|$\delta_i$|Inhibitory onset time(ms)|3-15|
|$\rho_i$|Inhibitory $tau$ ratio (fall/rise)|5-20|
|$\bar{g}_i$|Inhibitory max conductance|0.04-0.5|
|$\beta$|Inhibitory/Excitatory $tau$ rise ratio|0.5-5|

### Double exponential to explain the net synaptic conductance.


```python
alpha = exp(-(t-delta_e)/e_d) - exp(-(t-delta_e)/e_r)
```


```python
alpha
```


```python
#alpha = alpha.subs(e_d, (rho_e*e_r)).doit()
```


```python
alpha_prime = alpha.diff(t)
```


```python
alpha_prime
```


```python
theta_e = solve(alpha_prime,t) # Time to peak
```


```python
theta_e = logcombine(theta_e[0])
```


```python
theta_e
```


```python
N(theta_e.subs(averageEstimateDict))
```


```python
alpha_star = simplify(alpha.subs(t, theta_e).doit())
```


```python
alpha_star
```


```python
#alpha_star = simplify(alpha) # Replacing e_d/e_r with tau_e
```

### Finding maximum of the curve and substituting ratio of taus


```python
alpha_star
```


```python
E = Piecewise((0, t < delta_e), (g_e * (alpha/alpha_star), True))
```

### Final equation for Excitation normalized to be maximum at $g_e$


```python
E
```

### Verifying that E Behaves


```python
E_check = N(E.subs(averageEstimateDict))
```


```python
E_check.free_symbols
```


```python
plot(E_check,(t,0,100))
```

### Doing the same with inhibition


```python
I = E.xreplace({g_e: g_i, rho_e: rho_i, e_r:i_r, e_d: i_d, delta_e: delta_i})
```


```python
I
```

### Verifying that I Behaves


```python
I_check = N(I.subs(averageEstimateDict))
```


```python
plot(-I_check,(t,0,100))
```

### Now finding the control peak using difference of these double-exponentials


```python
C =  E*() - I
```


```python
C
```


```python
# 1. Get rid of the normalization factor.
# 2. The functions must approximate each other when close to 0.
# 3. The precision of the parameters in the equations must increase with g_e.
```

### Substituting excitatory and inhibitory ratios and putting $\delta_e$ to zero.


```python
#C = C.subs({g_i: g_e*P, i_r : e_r*b}) # Replacing g_i with P*ge
C = C.subs({delta_e:0})
```


```python
C_check = N(C.subs(averageEstimateDict))
```


```python
C_check
```

### Verifying that C behaves


```python
plot(C_check,(t,0,100))
```


```python
#C_check = N(C.subs({rho_e:7, rho_i: 15}))
```


```python
C_check
```


```python
C_check = C_check.subs(averageEstimateDict)
```


```python
plot(C_check,(t,0,100))
```


```python
C_prime = diff(C,t)
```


```python
C_prime_check = N(C_prime.subs(averageEstimateDict))
```


```python
plot(C_prime_check,(t,0,100))
```


```python
C_prime_prime = diff(C_prime,t)
```


```python
C_prime_prime_check = N(C_prime_prime.subs(averageEstimateDict))
```


```python
plot(C_prime_prime_check,(t,0,100))
```

<div></div>


```python
#simplify(C.subs(t,log(T)))
```


```python
#C.subs(delta_i, 1/g_e)
```


```python
#x, denominator = cse(simplify(C_prime.as_numer_denom()))
```


```python
#T = symbols('T')
```


```python
simplify(C_prime)
```

#### Explicit solving this equation doesn't work

### Trying to use lambert function


```python
a,b,c,d = -t/e_r, -t/e_d, -(t - delta_i)/i_r, -(t - delta_i)/i_d
```


```python
alpha_star
```


```python
W = lambda z: z*exp(z)
```


```python
LambertW(W(a))
```


```python
C = g_e*(exp(LambertW(W(a))) - exp(LambertW(W(b))))/alpha_star - g_i*(exp(LambertW(W(c))) - exp(LambertW(W(d))))/alpha_star.subs({ e_r:i_r, e_d: i_d})
```


```python
C
```


```python
C.diff(t)
```


```python
t_star = solve(expand(C.diff(t)),t)
```


```python
a
```


```python
exp(a).diff(t)
```


```python
mpmath.lambertw(1)
```




    mpf('0.56714329040978384')




```python
C.diff(t)
```


```python
powsimp(C.subs({-t/e_r:a, -t/e_d:b, -(t - delta_i)/i_r:c, -(t - delta_i)/i_d:d}))
```


```python

```


```python
piecewise_C_star = simplify(ratsimp(factor(C_prime))).args
```


```python
C_star_1 = simplify(piecewise_C_star[1][0])
```


```python
C_star_1.args
```


```python
simplify(solveset(C_star_1.args[5],t).doit())
```


```python
factor(C_star_1).collect(exp(t)).args
```


```python
expand_log(factor(C_prime))
```


```python
denominator
```


```python
-x[5]/((i_d*exp(x[13]*x[9])) - (i_r*(exp(x[13]/i_d))))
```


```python
j,k = symbols({'J','K'})
```


```python
new_eq = simplify(C_prime.subs({e_d:e_r*((j+1)/(j-1)), i_d:i_r*((k+1)/(k-1))}))
```


```python
refine(powsimp(new_eq.as_numer_denom()))
```


```python
eq_1 = latex("\bar{g}_e*\tau_{er}*(-P*\rho_e*\rho_i**(\rho_i/(\rho_i - 1) - 1 + 1/(\rho_i - 1))*(\rho_e**(\rho_e/(\rho_e - 1)) - \rho_e**(1/(\rho_e - 1)))*(\rho_i*exp((\delta_i - t)/(\beta*\tau_{er})) - exp((\delta_i - t)/(\beta*\rho_i*\tau_{er})))*exp(t*(\rho_e + 1)/(\rho_e*\tau_{er}))/(\beta*(\rho_i**(\rho_i/(\rho_i - 1)) - \rho_i**(1/(\rho_i - 1)))) + \rho_e**((\rho_e + 1)/(\rho_e - 1))*(\rho_e*exp(t/(\rho_e*\tau_{er})) - exp(t/\tau_{er})))")
```


```python
eq_1
```


```python
C.diff(t).diff(t)
```


```python
t_star = (delta_i - (b*e_r*log((P*(rho_i+1)*rho_e)/((rho_e+1)*rho_i))))/(b+1)
```


```python
t_star.subs(averageEstimateDict)
```


```python
N(delta_i.subs(averageEstimateDict))
```

<div></div>

### Unfortunately this is not possible: Since the $\tau_{decay}$ will not contribute to the first peak, we can eliminate them.


```python
erise = Piecewise((0, t < delta_e), (g_e * (exp(-(t-delta_e)/e_r)/alpha_star), True))
```


```python
efall = Piecewise((0, t < delta_e), (g_e * (exp(-(t-delta_e)/e_d)/alpha_star), True))
```


```python
irise = erise.subs({g_e: g_i, rho_e: rho_i, e_r:i_r, e_d: i_d, delta_e: delta_i})
```


```python
ifall = efall.subs({g_e: g_i, rho_e: rho_i, e_r:i_r, e_d: i_d, delta_e: delta_i})
```


```python
C = C.subs({g_i: g_e*P, i_r : e_r*b}) # Replacing g_i with P*ge
C = C.subs({delta_e:0})
```


```python
C
```


```python
C_check = C.subs({P:0, delta_i:2})
```


```python
C_check
```


```python
C
```


```python
C.diff(t)
```


```python
averageEstimateDict
```


```python
C_check = N(C.subs(averageEstimateDict))
```


```python
C_check
```

### Verifying that C behaves


```python
plot(C_check,(t,0,100))
```


```python
plot(erise.subs(averageEstimateDict),(t,0,100))
```


```python
plot(((efall-erise)-(-irise)).subs(averageEstimateDict),(t,0,100))
```


```python
plot(((efall-erise)-(ifall-irise)).subs(averageEstimateDict),(t,0,100))
```


```python
plot((E-I).subs(averageEstimateDict),(t,0,100))
```


```python
C_check = N(C.subs({rho_e:7, rho_i: 15}))
```


```python
C_check
```


```python
C_check = C_check.subs(averageEstimateDict)
```


```python
plot(C_check,(t,0,100))
```


```python
C_prime = diff(C,t)
```


```python
C_prime_check = N(C_prime.subs(averageEstimateDict))
```


```python
plot(C_prime_check,(t,0,100))
```


```python
C_prime_prime = diff(C_prime,t)
```


```python
C_prime_prime_check = N(C_prime_prime.subs(averageEstimateDict))
```


```python
plot(C_prime_prime_check,(t,0,100))
```

<div></div>


```python
simplify(C)
```


```python

```


```python

```


```python
C_prime = diff(C,t)
```


```python
C_prime
```


```python
theta_C = solve(C_prime, t)
```


```python
theta_C
```


```python
C_star = C.subs(t, theta_C[0])
```


```python
C_star = C_star.subs(delta_e,0.) # Putting excitatory delay to zero
```

### Assuming that certain ratios are more than one and substituting


```python
C_star = C_star.subs({i_d: (i_r*tau_i)+i_r, e_d: (e_r*tau_e)+e_r, g_i: g_e*P, i_r:e_r*b}) # Replacing e_d/e_r with tau_e
```


```python
C_star = cancel(powsimp(factor(C_star), deep=True))
```


```python
C_star = C_star.collect([g_e, delta_i, P])
```


```python
#C_star1 = limit(limit(C_star, (1/tau_e), 0), (1/tau_i),0)
```


```python
#simplify(C_star.subs({e_r: 4., i_r: 2.73, g_i:P*g_e, e_d: g_e*b, i_d : g_e*g, delta_i: theta/g_e }))
```


```python
#cancel(C_star1.together(deep=True))
```


```python
C_star.free_symbols
```


```python
#tau_e1, tau_i1 = symbols('\\tau_{e1} \\tau_{i1}', real=True, positive=True)
```


```python
#simplify(C_star.subs({tau_e:tau_e1+1, tau_i:tau_i1+1}))
```


```python
C_star = simplify(C_star)
```


```python
C_star
```


```python
cse(C_star)
```


```python
cse(simplify(diff(C_star,g_e)))
```


```python
x = Symbol('x')
```


```python
y = x**(1/x)
```


```python
y.subs(x,40).evalf()
```


```python
theta_C_nice = simplify(theta_C[0].subs({i_d: (i_r*tau_i)+i_r, e_d: (e_r*tau_e)+e_r, g_i: g_e*P}))
```


```python
cse(cancel(expand(theta_C_nice)))
```


```python
theta_C_nice
```


```python
limit(x/(x-1),x,5)
```


```python
log(-2)
```


```python

```

15a2914ebdf961755bf6d7b223dad29bd192fc41
ipynb
Jupyter Notebook
model/Model_double_exps.ipynb
elifesciences-publications/linearity
777769212ac43d854d23d5b967c6323747c56c09
MIT
2019-04-22T17:07:37.000Z
2019-04-22T17:07:37.000Z
model/Model_double_exps.ipynb
elifesciences-publications/linearity
777769212ac43d854d23d5b967c6323747c56c09
MIT
model/Model_double_exps.ipynb
elifesciences-publications/linearity
777769212ac43d854d23d5b967c6323747c56c09
MIT
2019-04-25T13:10:24.000Z
2021-09-05T03:45:36.000Z
Qwen/Qwen-72B
1. YES
2. YES