SharedSingle Neuron Equation / Single Neuron Equation.ipynbOpen in CoCalc
Ahhh

You should find the Single Neuron Equation assignment as a pdf on the course calendar at http://geofhagopian.net/m2c/M2C-S18/M2C-S18-Syllabus.pdf (under week 4 on the calendar). Use this space to start a collaboration with a partner on CoCalc, if you like. You can also use other means to complete the project.

Single Neuron Equation

Background theory

Biologists are currently applying mathematical modeling to a wide range of problems. One such problem involves
the function of nerve cells (neurons) in the brain. These cells generate electrical impulses that move along their
axons, affecting the activity of neighboring cells. Even for a single neuron, a complete model using what is currently
known about its firing mechanism may take several equations, Since the human brain conatains on the order of 101010^{10} neurons,
an exact model is not feasible. In this application you will be looking at a very simplified model for an isolated population
of neurons all having similar properties.

center>

Let x(t)x(t) denote the percent of neurons firing at time tt, normalized to be between 0 (low activity) and 1 (high activity).
A simple model representing the change of activity level in the population is given by the differential equation

where EE is the level of activity coming from the cells external to the population, θ\theta is a common threshold level for
cells in the set, and SaS_a is a response function that models the change in activity level due to a given input.

We will use a standard “sigmoidal” response function .

The nonlinear function SaS_a can be seen to increase monotonically from 0 to 1 as zz increases from to \infty. It
is called a sigmoidal function because it has a sort of stylized S-shape. You may remember that solutions of the logistic growth
equation had this same shape.

<image=neuron01.png>

  1. Find a formula for the derivative of Sa(z)S_a(z), and show that it satisfies the identity

*Hint: Mostly algebra.*
Completed on seperate sheet of paper

<image="neuron01.png">

  1. Draw a graph of Sa(z)S_a(z) for a=3,10,a=3,10, and 2020. Where is the slope a maximum? Is it the same in each case?

Explain how the graph of differs from the graph of Sa(z)S_a(z).
Because the value of SaS_a is always between 0 and 1, if x>1x>1 the slope is negative and
if x<0x<0 it is positive. This means that any equilibrium solutions, that is, values of xx where dxdt=0\dfrac{dx}{dt}=0
must lie between 0 and 1. It also implies that the arrows on the phase line are always pointing down above x=1x=1 and
up below x=0.x=0.

Figure 2.27 shows the graphs of y=xy=x (the dashed line) and the response function
for θ=0.4,0.7,\theta=0.4,0.7, and 1.01.0. It can be seen that for small θ\theta there will be one equilibrium solution near x=1x=1 and
for large θ\theta there will be one equilibrium near x=0x=0. This seems reasonable since a high threshold means it
takes a lot of input to produce very much activity. For θ\theta in the middle range, however, there can be three
equilibrium solutions.


*In Sage, use commands like
* `z,S=var('z S')`
`p0=plot()`
`p1=plot()`
`show(p0+p1)`
def findMax(a):
    locmaxz=0
    maxz=0
    z=0
z,S=var('z S')
a=3
t=0
p0=plot(1/(1+exp(-1*a*(z+t))))
a=10
p1=plot(1/(1+exp(-1*a*(z+t))))
a=20
p2=plot(1/(1+exp(-1*a*(z+t))))
show(p0+p1+p2)

(z-Theta) shifts the graph left or right

Slope is maximum at x = 0 for all values of a

3) Draw phase lines for with θ=0.4,0.5,,0.9,1.0\theta=0.4,0.5,\dots,0.9,1.0

Label each equilibrium point as a sink, source or node. You will need a numerical equation solver to find the
equilibrium values; that is, the value of xx where . As a check,
the three equilibria for θ=0.7\theta=0.7 are x10.007188,x2=0.5,x_1\approx 0.007188,x_2=0.5, and x30.992812.x_3\approx 0.992812.

*Hint To get some idea of where to hunt for what, it's helpful to have a graph of the seven curves for *
θ=0.4,0.5,,0.9,1.0\theta=0.4,0.5,\dots,0.9,1.0. In Sage,

g=Graphics()
x,y=var('x y')
p=plot(x,(x,-.2,1.2))
g=g+p
for i in range(4,11):
p=plot(1/(1+e^(-10*(x+0.2-i/10))),(x,-0.2,1.2))
g=g+p
show(g)
To find the numerical approximations in Sage, you can use the find_root command.

Hint: The easiest way to get the bifurcation curve is just to connect the dots by hand with a smooth curve.

g=Graphics()
x,y=var('x y')
for i in range(4,11):
    T = text("theta = "+str(i/10),(i/10-0.1,0.9-0.1*i))
    p=plot(-x+1/(1+exp(-10*(x+0.2-i/10))),(x,-0.2,1.2))
    g=g+p+T
show(g)
for i in range(4,11):
    try:
        h=find_root(-x+1/(1+exp(-10*(x+0.2-i/10))),0,0.1)
    except RuntimeError:
        continue
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y > 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y < 0:
            t = ("sink",(h))
            show(t)
        if y > 0:
            t = ("node",(h))
            show(t)
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y < 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y > 0:
            t = ("source",(h))
            show(t)
        if y < 0:
            t = ("node",(h))
            show(t)
for i in range(4,11):
    try:
        h=find_root(-x+1/(1+exp(-10*(x+0.2-i/10))),0.1,0.9)
    except RuntimeError:
        continue
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y > 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y < 0:
            t = ("sink",(h))
            show(t)
        if y > 0:
            t = ("node",(h))
            show(t)
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y < 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y > 0:
            t = ("source",(h))
            show(t)
        if y < 0:
            t = ("node",(h))
            show(t)
for i in range(4,11):
    try:
        h=find_root(-x+1/(1+exp(-10*(x+0.2-i/10))),0.8,1.2)
    except RuntimeError:
        continue
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y > 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y < 0:
            t = ("sink",(h))
            show(t)
        if y > 0:
            t = ("node",(h))
            show(t)
    f=h-0.1
    y=-f+1/(1+exp(-10*(f+0.2-i/10)))
    if y < 0:
        f=h+0.1
        y=-f+1/(1+exp(-10*(f+0.2-i/10)))
        if y > 0:
            t = ("source",(h))
            show(t)
        if y < 0:
            t = ("node",(h))
            show(t)
(sink,0.0224011111715)\left(\verb|sink|, 0.0224011111715\right)
(sink,0.00718806418267)\left(\verb|sink|, 0.00718806418267\right)
(sink,0.00253596891291)\left(\verb|sink|, 0.00253596891291\right)
(sink,0.000919458826073)\left(\verb|sink|, 0.000919458826073\right)
(sink,0.000336480036919)\left(\verb|sink|, 0.000336480036919\right)
(source,0.328504174745)\left(\verb|source|, 0.328504174745\right)
(source,0.5)\left(\verb|source|, 0.5\right)
(source,0.671495825255)\left(\verb|source|, 0.671495825255\right)
(sink,0.999663519963)\left(\verb|sink|, 0.999663519963\right)
(sink,0.999080541174)\left(\verb|sink|, 0.999080541174\right)
(sink,0.997464031087)\left(\verb|sink|, 0.997464031087\right)
(sink,0.992811935817)\left(\verb|sink|, 0.992811935817\right)
(sink,0.977598888829)\left(\verb|sink|, 0.977598888829\right)
  1. Make a bifurcation diagram by putting the seven phase lines from problem (3) in a row and joining the equilibrium points.
L1=line( [ (0.4, 0.999663), (0.5, 0.99908) ])
L2=line( [ (0.5, 0.99908), (0.6, 0.99746) ])
L3=line( [ (0.5, 0.99908), (0.6, 0.328504) ])
L4=line( [ (0.5, 0.99908), (0.6, 0.022401) ])
L5=line( [ (0.6, 0.99746), (0.7, 0.992817) ])
L6=line( [ (0.6, 0.328504), (0.7, 0.5) ])
L7=line( [ (0.6, 0.022401), (0.7, 0.007188) ])
L8=line( [ (0.7, 0.992817), (0.8, 0.97759) ])
L9=line( [ (0.7, 0.5), (0.8, 0.671495) ])
L10=line( [ (0.7, 0.007188), (0.8, 0.0025359) ])
L11=line( [ (0.8, 0.0025359), (0.9, 0.00033648) ])
L12=line( [ (0.8, 0.671495), (0.9, 0.00033648) ])
L13=line( [ (0.8, 0.97759), (0.9, 0.00033648) ])
L14=line( [ (0.9, 0.00033648), (1.0, 0.00091945) ])
show(L1+L2+L3+L4+L5+L6+L7+L8+L9+L10+L11+L12+L13+L14)
  1. From the bifurcation diagram, estimate the two bifurcation values of θ\theta where the number of equilibrium points
    changes from one to three, and then from three back to one.

*Hint: I think it’s better to use the figure showing the functions *
superimposed on y=xy=x

h=find_root(-1+(10*exp(-10*(x+0.2)))/((1+exp(-10*(x+0.2)))^2),0.0,1.0)
#we set the derivates equal to each other, and set the theta equal to 0
f=(1/(1+exp(-10*(h+0.2))))
#we then plug it into the origional
j=f-h
show(j)
#we then subtract the root from the new y value(j)
h=find_root(-1+(10*exp(-10*(x+0.2)))/((1+exp(-10*(x+0.2)))^2),-10,0)
#we set the derivates equal to each other, and set the theta equal to 0
f=(1/(1+exp(-10*(h+0.2))))
#we then plug it into the origional
j=f-h
show(j)
#we then subtract the root from the new y value(j)
0.8809546277311860.880954627731186
0.5190453722688140.519045372268814
These are the two bifurcation points
  1. Find the two bifurcation values of θ\theta analytically. You will need to solve two simultaneous equations obtained
    by using the fact that at a bifurcation value of θ\theta the curves and y=xy=x have a point
    of tangency. At this point the yy-values are equal and the slopes are also equal.

Hint: Teasing a numerical solution to a system of non-linear equations out of Sage is not straight-forward.
Wolfram Alpha chugs away at it for a while and then gives up. You may resort to guess and check...we already know about where

to look, after all. There is a clever way, I think.

See previous answer
  1. Use your computer algebra system to draw a slope field for the ode with a=10,E=0.2,a=10,E=0.2, and θ=0.7\theta=0.7.
    Let tt vary from 0 to 20. Use initial values x(0)=0.1,0.3,0.5,0.7,x(0)=0.1,0.3,0.5,0.7, and 0.90.9. Describe what happens to the
    activity as tt\rightarrow\infty.

Hint: In order to accomplish this in Sage, it's helpful to write the code for the building the numerical solutions.
The code below was adapted from the Mendel University (located in Brno, Czech Republic) site http://user.mendelu.cz/marik/sage/dr.pdf.
In Czeck, "krok" must mean something like "step size":

t,x=var('t x')
def rk(fnc,t0,x0,krok,t1):
g=fast_float(fnc,'t x','x')
n=int((1.0)*(t1-t0)/krok)
t00=t0; x00=x0
soln=[[t00,x00]]
for i in range(n+1):
m1=g(t00,x00)
m2=g(t00+krok/2,x00+m1*krok/2)
x00=x00+krok*m2
t00=t00+krok
soln.append([t00,x00])
return soln

(t_min,t_max,x_min,x_max)=(0,20,0,1)
f(t,x)=-x+1/(1+e^(-10*(x+0.2-0.7)))
A=plot_slope_field(f(t,x),(t, t_min, t_max), (x, x_min, x_max))
B=list_plot(rk(f(t,x), 0, 0.1, 0.1, t_max), plotjoined=True)
C=list_plot(rk(f(t,x), 0, 0.3, 0.1, t_max), plotjoined=True)
D=list_plot(rk(f(t,x), 0, 0.5, 0.1, t_max), plotjoined=True)
E=list_plot(rk(f(t,x), 0, 0.7, 0.1, t_max), plotjoined=True)
F=list_plot(rk(f(t,x), 0, 0.9, 0.1, t_max), plotjoined=True)
show(A+B+C+D+E+F,ymin=x_min,ymax=x_max,xmin=t_min,xmax=t_max)

This is somewhat shy of rk4, but it seems to work well enough and we get the slope field A on which we can use list_plot to superimpose the 5 solutions.

t,x=var('t x')
def rk(fnc,t0,x0,krok,t1):
    g=fast_float(fnc,'t x','x')
    n=int((1.0)*(t1-t0)/krok)
    t00=t0; x00=x0
    soln=[[t00,x00]]
    for i in range(n+1):
        m1=g(t00,x00)
        m2=g(t00+krok/2,x00+m1*krok/2)
        x00=x00+krok*m2
        t00=t00+krok
        soln.append([t00,x00])
    return soln

(t_min,t_max,x_min,x_max)=(0,20,0,1)
f(t,x)=-x+1/(1+e^(-10*(x+0.2-0.7)))
A=plot_slope_field(f(t,x),(t, t_min, t_max), (x, x_min, x_max))
B=list_plot(rk(f(t,x), 0, 0.1, 0.1, t_max), plotjoined=True)
C=list_plot(rk(f(t,x), 0, 0.3, 0.1, t_max), plotjoined=True)
D=list_plot(rk(f(t,x), 0, 0.5, 0.1, t_max), plotjoined=True)
E=list_plot(rk(f(t,x), 0, 0.7, 0.1, t_max), plotjoined=True)
F=list_plot(rk(f(t,x), 0, 0.9, 0.1, t_max), plotjoined=True)
show(A+B+C+D+E+F,ymin=x_min,ymax=x_max,xmin=t_min,xmax=t_max)
If t goes to infinity the f(t) reaches an equilibrium
  1. Redo problem (7) with periodic input E=E(t)=0.2(1+sin(t))E=E(t)=0.2(1+sin(t)). With this time varying input the equation
    is no longer autonomous. Explain carefully how the activity differs from that described in (7). Do you think
    there is a periodic solution separating the two types of solutions in the periodic case? This is an interesting
    problem and it might be a good time to look at the paper “Qualitative tools for studying periodic solutions
    and bifurcations as applied to the periodically harvested logistic equation” by Diego Benardete, V.W.Noonburg,
    and B. Pollina, Amer. Math. Monthy, vol 115, 202-219(2008). It discusses, in a very readable way, how one
    goes about determining the answer to such a question.

*HINT: At first it may seem simple enough to replace the 0.2 in the exponent of xx' with the sinusoid in Sage,
but as you may have come to expect, you get a cascade of red error messages. This caused me to look at the help file
for the fast_flot() function and modify it's call as shown below:

t,x=var('t x')
(t_min, t_max, x_min,x_max)=(0,20,0,1)
def rk(fnc,t0,x0,krok,t1):
g=fast_float(fnc,'t','x') #,'x')
n=int((1.0)*(t1-t0)/krok)
t00=t0; x00=x0
soln=[[t00,x00]]
for i in range(n+1):
m1=g(t00,x00)
m2=g(t00+krok/2,x00+m1*krok/2)
x00=x00+krok*m2
t00=t00+krok
soln.append([t00,x00])
return soln

f(t,x)=-x+1/(1+e^(-10*(x+0.2*(1+sin(t))-0.7)))
A=plot_slope_field(f(t,x),(t, t_min, t_max), (x, x_min, x_max))
B=list_plot(rk(f(t,x), 0, 0.1, 0.1, t_max), plotjoined=True)
C=list_plot(rk(f(t,x), 0, 0.3, 0.1, t_max), plotjoined=True)
D=list_plot(rk(f(t,x), 0, 0.5, 0.1, t_max), plotjoined=True)
E=list_plot(rk(f(t,x), 0, 0.7, 0.1, t_max), plotjoined=True)
F=list_plot(rk(f(t,x), 0, 0.9, 0.1, t_max), plotjoined=True)
show(A+B+C+D+E+F,ymin=x_min,ymax=x_max,xmin=t_min,xmax=t_max)

This should work!

t,x=var('t x')
(t_min, t_max, x_min,x_max)=(0,20,0,1)
def rk(fnc,t0,x0,krok,t1):
    g=fast_float(fnc,'t','x') #,'x')
    n=int((1.0)*(t1-t0)/krok)
    t00=t0; x00=x0
    soln=[[t00,x00]]
    for i in range(n+1):
        m1=g(t00,x00)
        m2=g(t00+krok/2,x00+m1*krok/2)
        x00=x00+krok*m2
        t00=t00+krok
        soln.append([t00,x00])
    return soln

f(t,x)=-x+1/(1+e^(-10*(x+0.2*(1+sin(t))-0.7)))
A=plot_slope_field(f(t,x),(t, t_min, t_max), (x, x_min, x_max))
B=list_plot(rk(f(t,x), 0, 0.1, 0.1, t_max), plotjoined=True)
C=list_plot(rk(f(t,x), 0, 0.3, 0.1, t_max), plotjoined=True)
D=list_plot(rk(f(t,x), 0, 0.5, 0.1, t_max), plotjoined=True)
E=list_plot(rk(f(t,x), 0, 0.7, 0.1, t_max), plotjoined=True)
F=list_plot(rk(f(t,x), 0, 0.9, 0.1, t_max), plotjoined=True)
show(A+B+C+D+E+F,ymin=x_min,ymax=x_max,xmin=t_min,xmax=t_max)
There are periodic solutions because the origional function has a sin function, which is periodic