Continuous models

This week, we first show how to go from a discrete-time to a continuous-time formulation of a model, namely, from difference to differential equations. We will get used to continuous-time models by building and analyzing a few of them. We will then learn how to integrate differential equations in a computer and eventually how to simulate a dynamics unfolding in continuous time.

In the next clip, LHD shows how translate a discrete-time SIS model in continuous time.

Derived the continuous-time model, he next shows how to find the equilibria and determine their stability.

SIS and SI dynamics are pretty special nonlinear models, for they are simple enough to be solved exactly (specifically, we can find a closed expression for I(t) valid at all times; see Bonus content below). In most cases – and "most" is an euphemism here – we are not able to do it, and our only possibility to watch the modeled system evolving is integrating the equations numerically. In the next clip, LHD introduces two basic methods of numerical integration (Euler's and Heun's methods). In class, we will see how to implement such methods in Python.

Whether we want to explore the behavior of a stochastic system or test how well our model does in predicting that behavior, we need a way to simulate the dynamics. While in discrete time we know when any of the next events may occur – we only have to test whether or not it does –, in continuous time we have a continuous distribution of times at which the next event could take place. (This distribution exists also in discrete time, but it has a trivial form: a single peak at time .) In the clip below, LHD shows how to sample such distribution.

Fast conversion from discrete to continuous

We have seen in the first clip how to go from a discrete-time description of the SIS model to its continuous-time counterpart. The procedure is the same for every model and in Sayama Ch. 6.3 you can find a general to connect the two. Here I want to show you an alternative way.

Consider the following generic equation for your discrete model,

where is the change in during the time step of length . Notice that depends on the length of the time step and we made such dependence explicit. The idea is to obtain the continuous formulation by expanding in a Taylor series around (assuming is differentiable, which generally is). We get

(at rigor, that derivative should be indicated as a partial derivative, i.e., ). First, , since no change can happen if no time has passed. Second, plugging Equation (2) back into Equation (1), substracting from both sides, and dividing everywhere by , we have

Taking the limit , we finally obtain

In summary, given the function quantifying the state change in your discrete model, take its derivative with respect to and then evaluate it at to obtain the corresponding continuous model.

Solving the SI(S) model

Let us briefly rederive the equation for the SIS model in the well-mixing or mean-field approximation – it is a good exercise. We assume that each of the individuals interacts with random individuals per time unit (e.g., per hour). Each of the contacts of each of the infected individuals in the population has probability of being susceptible. New infections are thus produced at rate , while infected individuals recover at rate . Considering a closed population (hence ), we get the following dynamic equation,

Before going further, as an exercise, try to derive Equation (4) starting from the respective discrete model and applying Equation (3). Given the more general definition of the model considered here, notice that the exponent of will not be , but , the expected number of infectious contacts (the model considered in the first clip is the special case ).

Back to Equation (4), rescaling via the mapping , we get the equivalent, but more coincise equation

Notice that an SI dynamics is obtained by setting (no recovery). The equation above is a second-order ODE of the family of Bernoulli differential equations. This kind of ODEs can be solved exactly by transforming them into linear ODEs. Here is the trick. Define , hence , and substitute in the equation above. Notice , from which,

Multuplying both sides by , we get the sought linear equation,

We have already seen how to solve this (reminder: solve first the homogeneous equation with no constant term, then solve the full equation by letting the factor in front of the exponential you got to depend on time); we get the solution

From , multiplying numerator and denominator by , we get the solution

First, notice that setting provides the solution for the SI model. Second, considered , we can recast the expression above in terms of the basic reproduction number, . Recall that is the equilibrium point the system reaches when (i.e., ). We can thus express the solution as

Over time, approaches if or vanishes if . In the former case, describes a sigmoid over time, with an initial exponential growth followed by an exponential saturation (see Sec. 2.2.2.2 here for proof); in the latter, just decays exponentially.

...

What if the system is exactly at the critical point (i.e., )? We cannot rely on the solution above, for we implicitly aasumed when solving Equation (6). Taking Equation (5) and setting , we still get a Bernoulli equation,

Proceeding as before, one obtains the solution

Therefore, the epidemic dies out also for , yet not exponentially, but hyperbolically () – at a much lower pace. This phenomenon of a (qualitatively!) longer relaxation time towards the equilibrium state is an example of critical slowing down. In particular, we see from Equation (7) that the time needed for the system to relax diverges linearly with the system's size – doubling requires doubling to reach the same value of .

Transversality of simple models

Can you think of other cases than infectious diseases where the dynamics is "isomorphic" to the SI(S) model? A couple of them come straight to my mind: a word or a behavior diffusing in a society, or a species growing in an environment. Wait...what?! How phenomena such different among them can be represented by the same model? Well, it turns out, they are not that different. Looking closer, they are all examples of population growth. The population of people using a word or adopting a behavior, or the population of a species. Words and behaviors "reproduce" themselves by being copied from mind to mind, through learning and imitation; species reproduce in the usual biological sense. This kind of transversality is typical of simple models: Stripped of finer details, many systems appear formally equivalent. You often want sophisticated models to make quantitatively accurate predictions, but simple, idealized models allow you to appreciate general patterns, build intuition and understanding, and transfer knowledge between fields.

After all these words, it's time to see a concrete example of such transversality. Consider a species in an environment with limited resources. Let be the abundance of the species at time and be the rate at which each individual reproduces. With no restriction on resource availability, the population would follow a never-ending exponential growth as the solution of the equation . Notice that, the same equation describes the more realistic situation where individuals are not ethernal beings and die at some rate . Just redefine via the mapping and you get the same equation above, with the difference that can now be negative too and so lead to an exponential extinction. This is Malthusian growth model. Another big step towards reality is to account for the fact that resources are never unlimited, and therefore individuals must compete for them. How do we encode such competition in the equation? A competition event needs at least two individuals to take place and causes the suppression of the loser (assuming there are always a winner and a loser), which either die or is unable to reproduce. The term accounting for competition should thus be negative and proportional to . We get to an equation of the form

with . This is the Verhulst's growth model. The quantity is called carrying capacity and you can easily check that it is the value the abundance eventually converges to when . Putting a limit on the resources available, prevents the population to grow indefinitely.

Now, look again at Equation (8). Any bell ringing in your brain? If not, go back by a few paragraphs. Can you spot a similar equation? Substitute with and Equation (8) is now identical to Equation (5). The reason is that, by parameter rescaling, both these equations can be expressed in the following form

known as logistic equation. Remember the name, because its discrete-time version we are going to see next week will probably surprise you!

You may ask at this point: what is the limited resource in the SIS model when used to represent contagions – spread of pathogens, words, behaviors? Well, susceptible individuals. These are finite in number and infected individuals "compete" to infect them – the more the people I infect, the less the people available for you to infect.

Reading between the (math) lines

The limited resources we made reference to do not show up in the equations. We just mentioned them to justify the presence of that quadratic term. They are hidden there and it is a good exercise to make them appear explicitly. To make the Verhulst model a bit more general, suppose we now have two species of abundances and , competing for the same, finite resource of abundance . Calling and the rates at which and individuals respectively consume the resource, and assuming for now a closed system, we get the following system of equations

If we think of the three abundances in terms of biomass (measured for instance as the mass of the present total organic carbon), we can directly compare them. Since the system is closed, the total biomass is constant over time (fast check: summing Equations (10) we get zero on the right hand side).

It seems like the equations for the two species do not talk to each other, as also the diagram seems to suggest. Where is the competition we stated at the beginning then? Let's not rush. Let us express in terms of and and define and . We get to the reduced system

Now read between the lines – Equations (11) have the same form of Equation (8), except for the extra term proportional to . It is the latter that accounts for the competition between the two species! Equations (11) are a specific instance of the so-called competitive Lotka-Volterra equations. We could have found them heuristically, as we did for Equation (8) (which, by the way, is the one-species version), instead we derived them from a more detailed model. As an exercise, make the system open by considering additional processes such as death or intake of new resources (or any other you can think of), and try to reduce the resulting system of equations. Do you still get terms of the form of those in Equations (11)? As a second exercise, find the condition under which Equations (11) reduce to two independent Malthusian equations.

Let us close this section by observing what does it mean in terms of closeness/openness of the system to reduce the number of equations. We already used the constraint to get rid of one equation (Equation (10c)). As a consequence, summing Equations (11) we no longer get zero on the right hand side, meaning the reduced system is not closed. We can see this graphically drawing the diagram corresponding to Equations (11).

The flow coming out from the and boxes, , is indeed the positive contribution to – lost biomass due to competition feeds back into the environment as new available resource –, while the rate in the loops, , give the negative contribution to it – species consume resources to reproduce.