RELATED ARTICLES

We therefore decided to develop a new course (LS30) from the ground up, that would take the approach of “modeling first.” Thus, we begin by asking the students to consider an ecosystem consisting of a predator species (“Sharks”) and a prey species (“Tuna”). We elicit from them what was likely to happen in such a system. (Fig. 1).

Fig. 1 On the first day of class, the students are presented with a simple ecosystem and are asked: what will happen to the Shark and Tuna populations?

Full size image

It quickly becomes clear that we need to make a model of this system to understand its long-term behavior.

Our focus is then on what mathematical concepts are necessary to make models of dynamical systems, particularly those that involve positive and negative feedback loops.

### How to Make a Model

We begin by developing the art of making a model, including learning how to choose state variables and how to represent their various feedforward and feedback relations to each other, both positive and negative. Then the state variables are collected into a geometric representation as a “state space” $${R}^{n}$$, in which the current state of a system is represented as a point in this state space, that is, a point in $${R}^{n}$$. (“…space is only a particular case of a triply extended magnitude.” (Riemann 1873) (Fig. 2).

Fig. 2 Geometry of “Shark-Tuna space.” Students learn in the first week that a point (T0, S0) represents the state of the system at a time

Full size image

In the case of the shark-tuna model, where the state variables are T (for the number of tuna) and S (for the number of sharks), the students then learn to ask: what makes T go up? What makes T go down? What makes S go up? What makes S go down? The various effects are summarized in a boxes-and-arrows flow diagram (Fig. 3 upper), and then definite mathematical terms are developed to make the effects precise. Time is spent on the art of modeling, and students learn why the tuna birth term is “bTT” and, most importantly, why the shark-meets-tuna term is “βST.” Then, when we develop models for chemical reactions, students are reinforced to think that, for example, “Sodium meets chloride” in a chemical reaction will be represented by “k[Na +][Cl-]” by exact analogy with predator and prey, or for that matter, with “Susceptible meets Infected” in an epidemiological model, where the key term is “βSI,” and β is the effectiveness of transmission, which reflects social distancing as well as inherent transmissibility.

Fig. 3 Upper: a typical model diagram for the shark-tuna problem. Lower: A “change equation” for the Shark-Tuna model. S’ and T’ are finite vectors

Full size image

The various effects, in their mathematical representations, are collected into what we call a “change equation” (Fig. 3 lower). Readers will recognize this as a multivariable first-order differential equation, but it is crucial to recognize that we have not yet defined the derivative. The left-hand side is the finite “change vector” (T’, S’) assigned to the point (T,S).

What we don’t do with a change equation is to treat it as an equation: we don’t try to use symbolic techniques for operating on its terms. We don’t try to substitute variables or use integration by parts, or any analytic techniques whatsoever. The change equation is there purely to set up a function from state space $${R}^{n}$$ to the space of all change vectors (“tangent space”), which is another copy of $${R}^{n}$$, but with different axes (S’ and T’ rather than S and T) (Fig. 4).

Fig. 4 Upper: the change equation is used to set up a vector field, a function from points in State Space to points in Tangent Space. Lower: At representative points in (T, S)-space, we draw the appropriate change vector (T’, S’)*1, where the constant 1 has the units of time

Full size image

In this way, we are replacing the vague and unhelpful nineteenth century concept of a “differential equation” with the twentieth century rigorous geometric concept of a vector field. The nineteenth-century concept, which is still the typical approach in traditional math Calculus courses, sees a differential equation as a piece of language, to be operated on by symbolic techniques in order to produce another piece of language called the “solution.” This is the style of the nineteenth century, which persists in today’s calculus classes. It is sobering to realize that if Cauchy rose from the dead this morning, he could be teaching freshman calculus this afternoon.

But differential equations in biology, even those as simple as the Shark-Tuna model, do not have analytical solutions. Therefore, it seemed pointless to us to develop paper-and-pencil techniques for “solving” differential equations, when those techniques would never be used in an actual application.

Instead, we emphasize the idea that state space is everywhere paved with change arrows, given by the change equation, and that the behavior of the system is generally obvious once those change arrows are visualized (Fig. 4 lower).

Students find this picture immediately appealing and intuitive: the state point follows the change arrows like a driver following directional instructions from Google Maps or Siri. As the students run their fingertips over the vector field, sketching out the trajectory, the process seems very natural, and it is useful to point out to them that what they are really doing is integrating a differential equation!

We found that the twentieth-century concept of vector field makes for better pedagogy than the nineteenth-century concept of differential equation. Vector fields were central to Poincare’s work “On the curves defined by a differential equation” (Poincaré 1892), which was the first work to look at differential equations geometrically. The ideas were heavily developed throughout the twentieth century, by pioneers like Birkhoff and Thom (See, for example, Abraham and Marsden’s Foundations of Mechanics (Abraham and Marsden 2008). The pedagogical value of the vector field concept was shown beautifully in Abraham’s picture book Dynamics-the Geometry of Behavior (Abraham and Shaw 1982).

Once we have a vector field, we can do the rigorous version of fingertip-tracing, and predict what behaviors will follow, by using Euler’s Method (Fig. 5 left) to approximate the true (unknown) solution curves numerically (Fig. 5 right).

Fig. 5 Left: Three steps of Euler’s method (blue arrows) to approximate the true (but unknowable) integral curve of the vector field (red curve). Right: The vector field for the Shark-Tuna model is shown in green. The blue trajectory was computed by many small steps of Euler’s method, starting from the initial condition given by the black dot

Full size image

### Programming

To employ Euler’s method, students learn a basic programming language in a weekly lab. We chose the Cocalc platform because it is free and because Python programming via Jupyter notebooks is readily available in this platform. Students perform a few steps of Euler’s method by hand and then write Python code to do many steps. Towards the end of the second quarter, they are writing Python code to calculate eigenvalues and eigenvectors of matrices. (Students report that “I know Python” is a powerful recommendation that helps them get into research labs, a goal for many of our undergrads.)

### Derivatives and Integrals

At this point in the course, we introduce the concept of the derivative. We treat one standard example where it is a “rate of change” and show informally and experimentally that the average rate of change around a point approaches a finite limit as the time interval shrinks. However, we do not then pursue technical lemmas, derivatives of famous functions or tricks for dealing with simple derivatives. Instead, we place great emphasis on using the concept of the derivative to define and compute the linear approximation to a function at a point. We do this because this interpretation is the true and important meaning of the derivative, one that is lost in most calculus classes, and also because it enables us to develop linear stability analysis of equilibrium points in 1D (Fig. 6). It also generalizes very naturally to the n-dimensional case: Toward the end of our 2-quarter sequence, we develop the notion of the Jacobian as the n-dimensional derivative, and then we can introduce linear stability analysis in n dimensions as the straightforward generalization of the 1D case.

Fig. 6 Linear stability analysis of the equilibrium points of the logistic vector field $${X}^{{\prime}}=bX(1-\frac{X}{k}$$). If the slope of the linear approximation is positive, the equilibrium point is unstable, while if it is negative, the point is stable

Full size image

We approach integrals as “the area under the curve” and then show using Riemann integration (“add up the little rectangles to get the area”) that the area under df(x)/dx is equal to f(x) and prove a simple version of the fundamental theorem of calculus. We focus on the concept of the integral, and not at all on special techniques for integrating elementary functions, integrals of famous functions, or any analytic procedures whatsoever for finding integrals of given functions. We are happy when our students can explain the relation between the COVID-19 “New Cases per day” graph and the “total cases” graph.

### Bifurcation

One of the most important contributions that dynamics can make to biology is to give us a rigorous foundation to the all-important notion of qualitative change. Qualitative changes are found throughout biology:

• cell fate “switches” (Ferrell and Machleder 1998)

• stem cell “commitment” (Park et al. 2012)

• ‘switches’ like the lysis/lysogeny switch in phage lambda (Ingalls, 2013, Garfinkel et al. 2017}

• “turning on” of the lac operon in e coli (Monod and Jacob 1961; Garfinkel et al. 2017)

• dormancy/germination “decision switch” in plants (FitzGerald and Keener 2021)

• ‘maturation’ of oocytes (Ferrell and Machleder 1998)

• Formation of “cell memory” (Lisman 1985)

• Emergence of cooperation between two competing species (Garfinkel et al. 2017).

• Emergence of an insect “outbreak,” for example, the Spruce Budworm (Ludwig et al. 1978; Garfinkel et al. 2017; Strogatz 2018)

The mechanism underlying these biological switches is a bistable dynamical system, generally arising out of a saddle-node bifurcation (Fig. 7).

Fig. 7 Left: A positive feedback loop “turns on” the lac operon: intracellular lactose activates production of lactose permease, which imports more lactose into the cell. Right: The positive feedback is shown in blue, while the red line represents metabolic degradation of lactose. Where the two curves cross are equilibrium points. Here the two outer equilibrium points are stable, and the middle one is unstable. This, then, is a classic bistable system: two stable equilibria flanking an unstable equilibrium (a “threshold”). As the parameter r changes, the system exhibits a saddle-node bifurcation

Full size image

The general phenomenon of switch-like behavior is not limited to biology and is also found in the social and psychological sciences. An interesting example is due to Robert Lux (Lux 1995: see also Garfinkel et al. 2017), who models the emergence of political polarization (a bistable system) in a previously centrist (monostable) system.

In his model, each person holds position Y (yes) or N (no). As our state variable, we will use X = the degree of imbalance in opinion, defined as the number of “Yes” opinions minus the number of “No” opinions. X will then change as people flip their opinions based on several factors, including a parameter a that reflects how sensitive the per capita flipping rate is to the degree of imbalance X. So the parameter a reflects the strength of the “bandwagon effect.”

The model undergoes a qualitative change as the parameter a changes: if a < 1, the system has a stable Equilibrium Point at X = 0 (balanced opinions). But if a > 1, the system undergoes a pitchfork bifurcation: the stable Equilibrium Point at X = 0 becomes unstable, and two new stable Equilibria appear at a positive and a negative value of X (political polarization) (Fig. 8).

Fig. 8 (upper) Dynamics of the Lux public opinion model for 3 different values of the bandwagon effect parameter a. For a < 1 (black), there is a single stable Equilibrium Point at X = 0, but for a > 1 (red and blue), the X = 0 EP becomes unstable, and two new stable EPs appear at positive and negative values of X. (lower) The location and stability of the EPs, plotted as a function of a, make a classic pitchfork bifurcation diagram

Full size image

### Oscillation

We then spend a full 2 weeks on the topic of oscillation in biology, and its description by limit cycle attractors in dynamical systems. We discuss many examples of oscillatory processes in biology, from oscillatory gene expression to insulin/glucose oscillations to gonadal hormone oscillations to body temperature rhythms (Fig. 9). The importance of these rhythmic processes in biology helps to emphasize the need for mathematics, because dynamical models are necessary to understand the mechanisms of oscillation in these systems.

Fig. 9 Oscillatory processes in physiology. Upper: oscillations in estradiol over 24 h, generated by a negative feedback loop in the Hypothalamus (Licinio et al. 1998). Middle: Oscillations in insulin and glucose over 20 h in a subject receiving constant IV nutrition (Sturis et al. 1991). Lower: oscillatory transcription of the gene p53 and its inhibitor Mdm2. (Lahav et al. 2004)

Full size image

In particular, we use dynamical models to show that the key factors causing oscillation in biological systems are steeply sloped negative feedback, acting in the presence of time delays (Fig. 10). The students learn to see this as a universal property of negative feedback systems, and they learn to see the deep analogy between shark/tuna dynamics and insulin/glucose dynamics (“Oh, it’s like the insulin eats the glucose”). We believe that this kind of analogical thinking is essential to the modeling process.

Fig. 10 Mechanism of oscillation in a negative feedback control system with a time delay due to an intermediate. The Hypothalamus–Pituitary–Gonadal hormone model exhibits oscillations for steeply sloped negative feedback (n = 10) but not for shallow slope feedback (n = 3)

Full size image

This completes the first quarter, comprising 10 weeks X 2.5 h/week.

### Chaos

In the second quarter, we spend a week on chaos. Students iterate the discrete logistic function (the May equation) and discover its properties in lab. They also study a 3-dimensional food chain model that exhibits chaos (Hastings and Powell 1991) (Fig. 11).

Fig. 11 Sensitive dependence on initial conditions in the Hastings/Powell food chain model. At t = 0, a tiny ball of 10,000 initial conditions (arrow) was placed in the “handle” region, and evolved forward in time by the model. Snapshots of subsequent states show the stretching and folding that is characteristic of chaotic attractors. By t = 2500, the 10,000 points have evolved to populate the attractor nearly uniformly

Full size image

### Linear Algebra

In the bulk of the second quarter, the course goes on to develop some key concepts of linear algebra: the concept of a linear function and its representation by a matrix. We do not develop many of the themes of a traditional linear algebra course, such as matrix inversion techniques, linear independence, subspaces, dimension, orthogonality, least-squares methods, etc. We take a minimal path to the concepts of eigenvalues and eigenvectors, which is our destination.

In keeping with the theme of the course, our real interest is in linear discrete time dynamical systems and their representation as iterated matrices. So we consider a linear transform represented by a matrix Mx, and we consider the discrete time series X, M(X), M(M(X)), …

It is easy to show the key lesson that the long-term behavior of an iterated matrix is domination by the principal eigenvector (Fig. 12).

Fig. 12 The long-term behavior of an iterated matrix dynamical system is domination by its principal eigenvector. A 2 variable linear dynamical system (Black Bear Juveniles and Adults) is given by an iterated matrix. Starting from a circle of initial conditions (black dots), one application of the matrix carries those points to the oval of points formed by the dark-grey dots, and a second application carries them to the light grey points. By the 5th iteration (red dots), the points are lying nearly on a straight line, outlining the principal eigenvector U

Full size image

We illustrate this point with an application the Google PageRank algorithm (Fig. 13). We develop the idea of a “points to” matrix, (aij = 1 if site j points to site i, else 0) and observe that what we want are sites (‘pages’) that are pointed to by sites that are pointed to by sites that are pointed to…. This is not an infinite regress, but a call for a repeated iteration of the “points to” matrix to discover its principal eigenvector, whose components are the weights of each page. The pages with larger weights in the principal eigenvector appear higher in the search engine results.

Fig. 13 A toy model “internet” with 4 “websites” gives rise to a 4 × 4 “points to” matrix. The principal eigenvector of this matrix gives the Google PageRank values

Full size image

### Partial Derivatives

Once the student knows what an n-dimensional linear function is, it is easy to define the n-dimensional derivative as the linear approximation to a multi-dimensional function at a point, generalizing the idea from our development of the derivative as a linear approximation in our discussion of single-variable calculus.

We show that the linear approximation to a multi-dimensional function at a point is given by the Jacobian matrix of partial derivatives (Fig. 14).

Fig. 14 A 2D surface Z = f(X, Y) is shown in green (left). At a point (X0, Y0), the linear approximation to f is shown as a blue tangent plane (right). A 2D vector field V(X, Y) = (f(X,Y), g(X,Y)) (where W = g(X, Y) is the second component of the vector field) then has its linear approximation at a point (X0, Y0) given by the Jacobian J

Full size image

We then do linear stability analysis in dimensions bigger than 1. In particular, we derive 2D Hopf bifurcations for a number of our negative feedback models, thereby demonstrating the mechanism of the oscillation (Fig. 15).

Fig. 15 A Hopf bifurcation in a genetic control model. The stable equilibrium point becomes unstable, and the system undergoes stable oscillations when the parameter k, reflecting the strength of the positive feedback loop, becomes greater than 10

Full size image