Function definitions and set up

In this exercise, we want to explore the dynamics of the Morris-Lecar model of the barnacle muscle fiber, which is similar to, but simpler than, the Hodgkin Huxley model of the squid giant axon. Here we are going use the theory of dynamical systems to examine the model geometrically, explore the phase portrait of the model and assessing the stability properties of the equilibria that represent asymptotic solutions to the model.

Recall that the Morris Lecar model involves two voltage gated channels, one each for calcium and potassium. If we fix the calcium gating variable \(m\) at its voltage dependent steady state \(m_{\infty}\), the model reduces to two differential equations, one for the voltage \(v\) and one for the potassium gating variable \(w\):

\[ C\frac{dv}{dt}=i-g_{Ca}m_\infty(v)(v-v_{Ca})-g_Kw(v-v_K)-g_L(v-v_L) \] and \[ \tau_w(v)\frac{dw}{dt}=\phi(w_\infty(v)-w) \] where \[ m_\infty(v)=0.5(1+\tanh(\frac{v-v_1}{v_2})) \] \[ w_\infty(v)=0.5(1+\tanh(\frac{v-v_3}{v_4})) \]

and \[ \tau_w(v)=\frac{1}{\cosh((v-v_3)/2v_4)} \]

To do this we are going to use functions that are defined in the file DMBpplane.r, and we are also going to use functions from the deSolve library.

library(deSolve)
source('/Volumes/kerkhoffa$/MathBioClass/DMBpplane.r') #REPLACE WITH YOUR OWN PATH to the file

Ellner and Guckenheimer present two sets of parameters (Table 5.1) that they use to illustrate the range of dynamical behvior that the model can exhibit. The two sets differ in the values of only two parameters, the calcium conductance \(g_{Ca}\) and the potassium the parameter \(\phi\), which determines the amplitude of the time constant for the potassium channel.

The Morris-Lecar model is already defined in the DMBpplane.r file, but we can redefine it here, just so it is in front of us. For the different functions we are using, we need two separate functions. the lowercase morrislecar is for vector field and nullcline plotting, while the uppercase MorrisLecar is for solving the equations and evaluating the Jacobian to assess stability. Note that we leave only the two parameters mentioned above to be supplied to the funtion, fixing all of the other parameters at the values specified by Ellner and Guckenheimer.

# morris-lecar model for phase arrow and nullcline plotting
# in this specification, we give initial values of v and w and a vector of the two parameters as
# input, and we get back a vector of the two equations
morrislecar=function(v,w,parms) {
  gca=parms[1]; phi=parms[2];  
  c((1/20)* (90-gca*0.5*(1+tanh((v+1.2)/18))*(v-120)-8*w*(v+84)-2*(v+60)),
    phi*(0.5*(1+tanh((v-2)/30))-w)*cosh((v-2)/60))
} 

# morris-lecar model for computing solution trajectories and evaluating stability
# in this specification, we give vectors of time, initial conditions, and parameters, and
# get back a list of the differential equations
MorrisLecar=function(t,y,parms) {
  gca=parms[1]; phi=parms[2];
  v=y[1]; w=y[2];
  dv=(1/20)* (90-gca*0.5*(1+tanh((v+1.2)/18))*(v-120)-8*w*(v+84)-2*(v+60));
  dw=phi*(0.5*(1+tanh((v-2)/30))-w)*cosh((v-2)/60);
  return(list(c(dv,dw)))
}

Exploring the Model: Lower conductance and amplitude

Now that we have defined the functions, we can explore the behavior of the two parameter sets. First let’s consider a situation with lower calcium channel conductance \(g_{Ca}=4.4\) and lower amplitude $=0.04. Since we have defined the functions, we can just supply the functions and appropriate parameters to the plotting functions defined in DMBpplane.r. The phasearrows function plots the vector field for the model, while the function nullclines draws the (yep, you guessed it!) nullclines, one for each variable.

plot(1,xlim=c(-60,40), ylim=c(0,0.6), type='n',xlab="v",ylab="w") #blank plot to add vectors and clines
phasearrows(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=25,parms=c(4.4,0.04), add=T)
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=250,parms=c(4.4,0.04), add=T)

Examine the plot above and think about the meaning of the vector field and the nullclines. The vector field shows the vectors tangent to the function at any point, giving us the trajectory for any combination of \(v\) and \(w\) in the phase space. The nullclines show combinations of \(v\) and \(w\) for which \(\frac{dv}{dt}\)=0 (in blue) and \(\frac{dw}{dt}\)=0 (in green). Notice how the vectors tend towards horizontal around the \(v\) nullcline and towards vertical around the \(w\) nullcline. Why does this make sense? Where the nullclines cross, we find the equilibrium for this model.

To examine the properties of the equilibrium, we can use the newton function to find the equilibrium and compute its Jacobian matrix. The newton function also plots the point in the phase space, using symbols that describe its properties. Buried in the newton code is the following explanation: “Plotting symbol is closed/open=stable/unstable, circle/triangle=eigenvalues imaginary/real.” So let’s see what we have. We will first plot the nullclines again, then look at the results of the analysis.

plot(1,xlim=c(-60,40), ylim=c(0,0.6), type='n',xlab="v",ylab="w")
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=250,parms=c(4.4,0.04), add=T)
newton(MorrisLecar,x0=c(-30,0.1),parms=c(4.4,0.04))

## 1 -26.53938 0.1275774 
## 2 -26.59739 0.1293757 
## 3 -26.59687 0.1293793 
## Fixed point (x,y) =  -26.59687 0.1293793 
## Jacobian Df= 
##              [,1]         [,2]
## [1,] 0.0258199719 -22.96125321
## [2,] 0.0003351416  -0.04462988
## Eigenvalues 
## [1] -0.00940496+0.08033975i -0.00940496-0.08033975i
## $x
## [1] -26.5968670   0.1293793
## 
## $df
##              [,1]         [,2]
## [1,] 0.0258199719 -22.96125321
## [2,] 0.0003351416  -0.04462988

This analysis confirms that we have a single equilibrium under this set of parameters, located at (-26.6, 0.129), where the nullclines intersect. Since both eigenvalues are negative, we have a stable equilibrium, so it should attract trajectories that are in its neighborhood. However, we could still observe periodic orbits, so we probably want to explore some trajectories starting close to the equilibrium and maybe a bit further away. To do this, we can set the initial conditions and solve the model using the funciont lsoda. By playing around with different initial conditions, we can find a point where a small change in our starting value switches us from a trajectory that converges to the fixed point to a large periodic orbit around it.

# solve the Morris Lecar model, given initial conditions using lsoda
times<-seq(0,300,0.2)
outMLaIn<-lsoda(y=c(-26,0.1135),times=times, MorrisLecar, parms=c(4.4,0.04))
outMLaPer<-lsoda(y=c(-26,0.1134),times=times, MorrisLecar, parms=c(4.4,0.04))

#plot in phase space x=voltage, y=w
plot(1,xlim=c(-60,40), ylim=c(0,0.6), type='n',xlab="v",ylab="w")
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=250,parms=c(4.4,0.04), add=T)
lines(outMLaPer[,2], outMLaPer[,3], xlim=c(-60,40),ylim=c(0.1,0.6), type="l", col="black", add=T)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
lines(outMLaIn[,2], outMLaIn[,3], xlim=c(-60,40),ylim=c(0.1,0.6), type="l", col="red", add=T)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter

The red and black curves in this plot provide a phase portrait of the model, tracing the trajectory of the system over time as it moves through the the phase space as \(v\) and \(w\) dynamically change. In the red trajectory, we can see that the system spirals in towards the equilibrium point. However, we we start with just slightly different initial conditions (\(w=0.1134\) instead of \(w=0.1135\)) we get a very different asymptotic behavior, with the system falling into a stable periodic orbit around the equilibrium. To see what that looks like over time, we can plot \(v\) and \(w\) as a function of time.

plot(outMLaPer[,1],outMLaPer[,2], type="l", xlab="time", ylab="v", col="blue")
par(new=TRUE)
plot(outMLaPer[,1],outMLaPer[,3], xlab="", ylab="", type="l", axes=FALSE, col="forestgreen")
axis(side=4)
mtext(side=4, line=2.5, 'w')

From this plot, we can see that the two variables are periodically fluctuating out of phase with one another. These fluctuations reflect the trace of the periodic orbit in the phase space. Note that the scales on the variables are very different, which makes sense since the \(v\) is a voltage and \(w\) is a gating variable, i.e., a probability.

Biological question: The Morris Lecar model is specifically for the muscle fiber under a voltage clamp, rather than physiologically realistic conditions. Why build a mathematical model of membrane function for artificial, lab conditions?

Exploring the Model Further: Higher conductance and amplitude

If we change the underlying parameters of the model, the dynamics of the system can sometimes change substantially. Here, we raise the calcium conductance to \(g_{Ca}=5.5\) and the amplitude of the time constant to \(\phi=0.22\). We can again start by looking at the vector field and nullclines.

plot(1,xlim=c(-60,40), ylim=c(0,0.6), type='n',xlab="v",ylab="w") #blank plot to add vectors and clines
phasearrows(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=25,parms=c(5.5,0.22), add=T)
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=250,parms=c(5.5,0.22), add=T)

While the overall shape of the nullclines is similar, they now cross at three different points, so we need to examine the properties of all three equilibria. To find each equilibrium with Newton’s method, we have to start in the neighborhood of the equilibrium. Let’s begin with the leftmost, and move to the right

plot(1,xlim=c(-60,40), ylim=c(0,0.6), type='n',xlab="v",ylab="w")
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.6),resol=250,parms=c(5.5,0.22), add=T)
newton(MorrisLecar,x0=c(-30,0.1),parms=c(5.5,0.22)) #leftmost equilibrium
## 1 -22.70985 0.1515668 
## 2 -21.14534 0.1754925 
## 3 -21.09309 0.1766014 
## 4 -21.09315 0.1766017 
## Fixed point (x,y) =  -21.09315 0.1766017 
## Jacobian Df= 
##            [,1]        [,2]
## [1,] 0.18612985 -25.1627403
## [2,] 0.00229266  -0.2364972
## Eigenvalues 
## [1] -0.0251837+0.1141761i -0.0251837-0.1141761i
## $x
## [1] -21.0931494   0.1766017
## 
## $df
##            [,1]        [,2]
## [1,] 0.18612985 -25.1627403
## [2,] 0.00229266  -0.2364972
newton(MorrisLecar,x0=c(-10,0.28),parms=c(5.5,0.22)) #middle
## 1 -11.1487 0.2937578 
## 2 -11.50709 0.288929 
## 3 -11.51713 0.2888158 
## 4 -11.51714 0.2888157 
## Fixed point (x,y) =  -11.51714 0.2888157 
## Jacobian Df= 
##             [,1]        [,2]
## [1,] 0.453555938 -28.9931425
## [2,] 0.003089324  -0.2256066
## Eigenvalues 
## [1]  0.27443099 -0.04648161
## $x
## [1] -11.5171438   0.2888157
## 
## $df
##             [,1]        [,2]
## [1,] 0.453555938 -28.9931425
## [2,] 0.003089324  -0.2256066
newton(MorrisLecar,x0=c(4,0.5),parms=c(5.5,0.22)) #rightmost

## 1 3.028265 0.5171423 
## 2 2.975248 0.5162485 
## 3 2.97492 0.5162429 
## Fixed point (x,y) =  2.97492 0.5162429 
## Jacobian Df= 
##             [,1]       [,2]
## [1,] 0.372195448 -34.789968
## [2,] 0.003663281  -0.220029
## Eigenvalues 
## [1] 0.0760832+0.1994065i 0.0760832-0.1994065i
## $x
## [1] 2.9749200 0.5162429
## 
## $df
##             [,1]       [,2]
## [1,] 0.372195448 -34.789968
## [2,] 0.003663281  -0.220029

From the plotting symbols, it is immediately apparent that the three equilibria have different properties, and we can use the analysis of the Jacobian to figure out what they are. The first equilibrium (solid circle, -21.1, 0.177) has two negative eigenvalues and so is again a stable focus; it should attract trajectories in its neighborhood, as in the previous example.

The second equilibrium (open triangle, -11.5, 0.289) has one positive and one negative eigenvalue. We term this a saddle because trajectories will move toward it from some directions, but will tend to veer away, like a marble rollin on a horse’s saddle.

The third equilibrium (open circle, 2.97, 0.516) has two positive eigenvalues and is thus an unstable node. Trajectories in the neighborhood will tend away from this point.

We can see the varied nature of the equilibria if we start trajectories in the neighborhood of each.It is especially interesting to start very near, but on opposite sides of the saddle equilibrium.

# solve the Morris Lecar model, given initial conditions using lsoda
times<-seq(0,300,0.2)
outMLa1<-lsoda(y=c(-12,0.3),times=times, MorrisLecar, parms=c(5.5,0.22))
outMLa2<-lsoda(y=c(-10,0.3),times=times, MorrisLecar, parms=c(5.5,0.22))

#plot in phase space x=voltage, y=w
plot(1,xlim=c(-60,40), ylim=c(0,0.7), type='n',xlab="v",ylab="w")
nullclines(fun=morrislecar,xlim=c(-60,40),ylim=c(0,0.7),resol=250,parms=c(5.5,0.22), add=T)
lines(outMLa1[,2], outMLa1[,3], xlim=c(-60,40),ylim=c(0.1,0.6), type="l", col="black", add=T)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter
lines(outMLa2[,2], outMLa2[,3], xlim=c(-60,40),ylim=c(0.1,0.6), type="l", col="red", add=T)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a
## graphical parameter

To think about the differences between these two different trajectories, we can also plot \(w\) as a function of time for the two different initial conditions, to see how the behavior changes depending on only a 2 mV difference in the initial conditions.

plot(outMLa1[,1],outMLa1[,3], ylim=c(0,0.8), type="l", xlab="time", ylab="w", col="black")
lines(outMLa1[,1],outMLa2[,3], col="red")

Here we can see that just a 2 mV difference in initial conditions can greatly affect the asymptotic behavior of the model, moving it from dampening to strongly periodic oscillations. Remember that other than the initial conditions, everything else is exactly the same in this model.

Biological Question: What does this sort of sensitivity in simple models of biological function suggest about the range of behaviors possible for real biological phenomena?