Math 1241: Calculus and Dynamical Systems in Biology (Fall 2013)

Characterization of linear 2D flows

A tale of two eigenvalues

Any linear homogenous 2D system of ordinary differential equations may be written in the form $\left\{ \begin{array}{r c l} x' & = & \alpha x + \beta y \\ y' & = & \gamma x + \delta y \\ \end{array} \right.$ Setting $$x'=y'=0$$ we can easily find that the unique equilibrium point for a homogenous linear system is the origin, i.e. $$E=(0,0)$$ for the above system. In the next few sections we will characterize the different types of solutions for various parameters in order to later understand conditions for the orgin to be either a stable equilibrium point or an unstable equilibrium point.

growth and decay

In 2 dimensions the simplest linear system is really just 2 one dimensional equations which exhibit exponential growth or decay in each variable. In other words we have the following decoupled system $\left\{ \begin{array}{r c l} x' & = & ax \\ y' & = & by \;\;. \\ \end{array} \right.$ for the variables $$x=x(t)$$ and $$y=y(t)$$ with parameters $$a$$ and $$b$$. Since each of these equations is "decoupled" from the other we can solve them exactly just as we would for a 1D problem so that $x(t) = x(0)e^{at}$ $y(t) = y(0)e^{bt}.$

The Geogebra applet below plots this solution for different values of $$a$$, $$b$$ and initial condition $$\mathbf{X_0} = (x(0),y(0))$$.

If $$a>0$$ then we see exponential growth along the x-direction and if $$a<0$$ we see exponential decay along the x-direction and similarly the sign of the parameter $$b$$ controls growth and decay in the y-direction. This allows us to characterize three cases

(case I) exponential growth in both directions

(0,0) unstable equilibrium point

$$\lambda = +r_1, +r_2$$

(case II) exponential growth in one direction and exponential decay in the other direction

(0,0) unstable equilibrium point

$$\lambda = +r_1, -r_2$$

(case III) exponential decay in both directions

(0,0) stable equilibrium point

$$\lambda = -r_1, -r_2$$

In other words both $$a$$ and $$b$$ must be negative for the solution to converge to the equilibrium $$(0,0)$$.

pure oscillation

The system $\left\{ \begin{array}{r c l} x' & = & -a^2y \\ y' & = & b^2x \;\;. \\ \end{array} \right.$ exhibits elliptical orbits about the equilibrium $$E= (0,0)$$. The paramters above are squared to ensure that the sign of each coefficient does not change as the parameters $$a$$ and $$b$$ vary. In the geogebra applet below, we again start the dynamics at a location $$\mathbf{X_0} = (x(0),y(0))$$, but in this case will follow periodic orbits about the equilibrium. The green vector $$V$$ is the tangent vector to the curve at a given point.

We classify these elliptic orbits as stable. Even though the solution does not converge to the equilibrium point the solution "stays close."

oscillatory growth/decay

$\left\{ \begin{array}{r c l} x' & = & -y \\ y' & = & x + (2e)y \;\;. \\ \end{array} \right.$

We can characterize the behaviour of this system as follows:

(case IV) Oscillation and decay (Spiral in)

$$E= (0,0)$$ stable for $$e<0$$

$$\lambda = -\tau \pm i\omega$$

(case V) Oscillation and growth (Spiral out)

$$E= (0,0)$$ unstable for $$e>0$$

Simplified general linear system

The system $\left\{ \begin{array}{r c l} x' & = & -y \\ y' & = & -dx + ty \\ \end{array} \right.$ summarizes all the previous cases together. The two parameters $$d$$ and $$t$$ are sufficient to create an equivalent system to most any case for the original equations $\left\{ \begin{array}{r c l} x' & = & \alpha x + \beta y \\ y' & = & \gamma x + \delta y \\ \end{array} \right.$ with parameters $$\alpha$$, $$\beta$$, $$\gamma$$ and $$\delta$$. The applet below plots solutions for the two parameters $$d$$ and $$t$$ similar to the previous applets with the addition of plots of the Nullclines. The blue line is the x-Nullcline $$y=0$$ and the red line is the y-Nullcline $$d\cdot x=t\cdot y$$.

Can you determine conditions for $$d$$ and $$t$$ for stability of the equilibrium $$E= (0,0)$$?

Linear Algebra

In order to classify the behaviour of a linear flow in the simplest way we will use two ideas from Linear Algebra for a matrix. Any 2x2 matrix, $$A$$ can be written as $A = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \;\;.$ We will need to know how to find the determinant and the trace of such a matrix, both of which return a single number for any given matrix.

Trace

The trace of a matrix is simply the sum of the diagonal elements. For the 2x2 matrix $$A$$ above we have $tr(A) = a + d .$

Determinant

The determinant is a little more complicated than the trace. For the 2x2 matrix $$A$$ above we write its determinant as $\det(A) = |A| = ad-bc .$ The 2x2 determinant is simply the product of the diagonals minus the product of the off-diagonals.

Examples

Find the determinant and trace of the matrices P,Q,R below. \begin{align} P &= \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} & Q &= \begin{bmatrix} 3 & 0 \\ 1 & -2 \end{bmatrix} & R = \begin{bmatrix} 1 & -1 \\ 3 & 2 \\ \end{bmatrix} \\ \det(P) &= (1)(4)-(2)(3) & \det(Q) &= (3)(-2)-(0)(1) & \det(R) &= (1)(2)-(-1)(3) \\ \det(P) &= -2 & \det(Q) &= -6 & \det(P) &= 5 \\ tr(P) &= 1 + 4 & tr(Q) &= 3 -2 & tr(R) &= 1 + 2 \\ tr(P) &= 5 & tr(Q) &= 1 & tr(P) &= 3 \end{align}

Determining Stability of Nonlinear Equilibrium Point

When we studied the stability for a 1D ODE we found that the linearization determines the stability. Namely, if we have the dynamical system $x' = f(x)$ with $$f(x_e)=0$$ then

 $$f'(x_e) < 0$$ The equilibrium $$x=x_e$$ is stable $$f'(x_e) > 0$$ The equilibrium $$x=x_e$$ is unstable.

Using the previous two sections we outline a similar method to determine stability for the system $\left\{ \begin{array}{r c l} x' & = & f(x,y) \\ y' & = & g(x,y) \;\;. \\ \end{array} \right.$

Step 1 determine equilibria

First, we need to determine the equilibrium point of the system for which we would like to determine the stability. This can be done by solving $\left\{ \begin{array}{r c l} f(x_0,y_0) & = & 0 \\ g(x_0,y_0) & = & 0\\ \end{array} \right.$ for an equilibrium point $$(x_0,y_0)$$.

Step 2 compute the Jacobian

In the 2D case the derivative at a point is no longer just a number, as in the 1D case, but rather a 2x2 matrix. This is the matrix of derivatives, also called the Jacobian, that we calculated earlier in the semester. $J(x,y) = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix}$ Just like in the 1D case we need to evaluate the derivative at the equilibrium point. So let's set $A = J(x_0,y_0) = \begin{bmatrix} \frac{\partial f}{\partial x}(x_0,y_0) & \frac{\partial f}{\partial y}(x_0,y_0) \\ \frac{\partial g}{\partial x}(x_0,y_0) & \frac{\partial g}{\partial y}(x_0,y_0) \end{bmatrix}$

Step 3 calculate det(A) & tr(A)

We can determine the stabilty by calculating the determinant and trace of the derivative of the vector field evaluated at the equilibrium point. The last applet above plots solutions for the variable paramters $$t$$ and $$d$$ which are precisely the trace and determinant respectively. If the determinant is negative we are guaranteed to have exponential growth in one direction and exponential decay in another direction (not necessarily the x or y direction but an arbitrary direction in the plane.) So if the determinant is negative we have an unstable equilibrium point as in the case (II) under growth/decay. If the determinant is positive and the trace is positive we have either case III or case V, exponential growth with or without oscillations. If the determinant is positive and the trace is negative we have either case I or case IV, exponential decay with or without oscillations.

tr(A)$$<0$$ and det(A)$$>0$$

The equilibrium $$(x,y) = (x_0,y_0)$$ is stable. cases: I, IV

$$\Re(\lambda) < 0$$

tr(A)$$>0$$ or det(A)$$<0$$

The equilibrium $$(x,y) = (x_0,y_0)$$ is unstable. cases: II, III, V

$$\Re(\lambda) > 0$$

The case tr(A)=0, det(A)>0 is pure oscillation for the linear map. One may expect the nonlinear equilibrium to also exhibit pure oscillation but the first derivative is not enough to determine this sowe neglect this case in this class. tr(A)=0 is similar to the result f'(x_e)=0 in 1D. The case det(A)=0 is a degenerate case in which the 2D system was actually just a 1D system in disuise. We will not consider either of these cases in this class.

Need example...

Math 1241: Calculus and Dynamical Systems in Biology (Fall 2012)

The course webpage is maintained by the lecturer, Professor Duane Nykamp, here.

Catalog Description:
Differential and integral calculus with biological applications. Discrete and continuous dynamical systems. Models from fields such as ecology and evolution, epidemiology, physiology, genetic networks, neuroscience, and biochemistry.

Throughout the semester I will update the following outline of topics as we progress through the course. This will be a collection of topics which I feel are important. Something which is less important but perhaps of general interest is the cascades of the discrete logistic equation, explored below.

Period doubling of the logistic map

In class we looked at stability of the equilibria, $$E=0,\frac{1-a}{b}$$ for a discrete dynamical system with quadratic updating function \begin{aligned} %\left \{ x_{n+1} & = ax_n + bx_n^2 \\ x_0 & = x_0 \;\; . \end{aligned} Taking $$b=-a$$ this turns into the well known logistic map $$f(x_n)=ax_n(1-x_n)$$.

We don't consider parameter values in the range $$0< a <1$$ for this equation since the equilibrium $$\frac{1-a}{b}=\frac{a-1}{a}$$ is negative. For $$1< a <3$$ $$E=0$$ is unstable and $$E=1$$ is stable so any small positive initial population increases until it reaches the equilibrium $$E=1$$. This tells us that addition of the nonlinear term $$-ax_n^2$$ stablizes the population for growth rates in the regime $$1< a <3$$. But what happens for $$a >3$$?

Well for $$3< a <4$$ something interesting happens. For $$a >3$$ both equilibria are unstable, but this does not mean that solutions grow without bound. The image above is a plot of asymptotically stable periodic orbits . The image was generated numerically in Matlab using a simple code which calculates the first 1100 iterations of an orbit for a given parameter value and then plots the last 100 points of that orbit. This way we see where the solution is going to end up after a little time. As it turns out for $$3< a <1 + \sqrt{6}$$ the solution just bounces between two points. For example, if $$a=3.2$$ then, according to this figure, the orbit would be close to {0.8,0.5,0.8,0.5,0.8,...}. Instead of the solution geting closer and closer to a constant the solution gets closer and closer to oscillating between two constants for $$3< a <1 + \sqrt{6}\approx 3.45$$. For $$a$$ just above $$1 + \sqrt{6}$$ the period doubles so that after enough time the solution looks like {x_1,x_2,x_3,x_4,x_1,x_2,x_3,x_4,...} as the paramter $$a$$ is increased the period doubles again and again until about $$a=3.57$$. Another interesting observation is that after this doubling there are windows where small periods appear again, for example the largest window is just after $$a=3.8$$ where the system approaches a period 3 orbit.

We could also ask what happens for $$a>4$$ but if this is the case system will no longer be well defined. In the background we have assumed that $$x_n$$ is always in the interval $$[0,1]$$. Using a little calculus it's clear that $$x(1-x)$$ has a maximum of $$\frac{1}{4}$$ so that if $$a>4$$ there will be $$x_n$$ such that $$x_{n+1}>1$$ which causes the model to break down.