As in the case of the TBT module, the SIR module implements a set of interfaces that define the underlying process and advertise the ways to visualize the dynamics, see Fig. 14.1.

The capabilities of the SIR module are advertised and specified through the following interfaces:

Discrete

for discrete individual traits (or types, or compartments),

HasIBS

for individual based simulations, IBS,

HasDE.ODE

for deterministic dynamics based on ordinary differential equations, ODE,

HasDE.SDE

for stochastic dynamics based on the Euler-Maruyama method (ODE with a stochastic noise term),

HasDE.PDE

for deterministic dynamics in spatial settings modelled by partial differential equations, PDE, based on the Euler method combined with a reaction-diffusion term.

Moreover, the SIR module implements the following interfaces to request a set of GUI views to graphically visualize various aspects of the dynamics:

HasPop2D.Traits

display the types of all individuals in a 2D plot (IBS and PDE only),

HasPop3D.Traits

display the types of all individuals in a 3D plot (IBS and PDE only),

HasMean.Traits

plot the mean frequency of each type over time,

HasHistogram.Degree

plot the degree distribution of the population structure (IBS and PDE only),

HasHistogram.StatisticsStationary

plot the histogram for the sojourn time for the different type counts in the stationary distribution.

This tutorial guides you through the process of creating a new module from scratch but will eventually be identical to the existing SIR module. Since module names (as well as java class names) must be unique, we have to choose a name other than SIR. Let’s use SIRE – appending an E for EvoLudo.

The initial steps to create a new module are essentially the same as in the previous tutorial, chapter 13: Step-by-step tutorial $𝑰$: Re-creating the TBT module. In order to get started, create a new java file org.evoludo.simulator.modules.SIRE in the EvoLudoCore maven module. In VS Code right-click on the modules folder in EvoLudoCore/src/main/java/org/evoludo/simulator and select New Java File ¿ Class…. Enter the name of the new file as SIRE. The new file should look like this: The initial steps to create a new module are essentially the same as in the previous tutorial, chapter 13: Step-by-step tutorial $𝑰$: Re-creating the TBT module. In order to get started, create a new java file org.evoludo.simulator.modules.SIRE in the EvoLudoCore maven module. In VS Code right-click on the modules folder in EvoLudoCore/src/main/java/org/evoludo/simulator and select New Java File ¿ Class…. Enter the name of the new file as SIRE. The new file should look like this:

The following code snippets contain little to no comments because the SIR class serves as a reference file. Besides, the motivation and purpose of each code snippet is explained in detail in the accompanying text.

The new SIRE class must again extend the Discrete class because we are dealing with a discrete set of traits or types, say $S$, $I$, and $R$. This is done by changing the class definition from SIRE to SIRE extends Discrete. The Discrete class provides a number of useful methods and fields to handle discrete traits. The Discrete class is located in the same package as the SIRE class, so no import statement is needed.

However, extending Discrete triggers two errors in VS Code because (i) the Discrete class does not provide a default constructor and (ii) SIRE does not implement the abstract method getTitle() of the Discrete class. Both errors are easily addressed in VS Code by hovering over the error location and selecting Quick fix…, followed by Add constructor SIRE(EvoLudo) as well as Add unimplemented methods.

Now the content of the file looks as follows:

In the constructor, there is nothing else we need to do, so simply delete the TODO comment.

In essence, the $SIR$ model describes a contact process (durrett:XXX) with multiple types. In contrast to the $TBT$ tutorial, payoffs do not matter, and no game theoretical interactions take place between individuals. Instead, infected individuals may pass the infection to susceptible ones – should they come into contact. The other types can spontaneously convert to another type: infected ones recover, and recovered ones may lose resistance and become susceptible again. Since interactions with other types do not involve payoffs, SIRE does not implement the interface Payoffs and no further methods need to be implemented.

Let us start with some housekeeping. For convenience and to keep the code more readable, we define three constants for the indices of the different types:

Next, when loading a new module, the method load() is called and provides an opportunity to initialize the module. Most importantly, we must specify the number of types (or traits) by setting the field nTraits. Optionally, we may also set the names and colours of each trait. All this is done by overriding the method load():

The Color class requires another import statement at the top:

The classical $SIR$ model in the form of an ODE has already been derived in the introduction, see Eq. (14.2). So, let’s start with the implementation of the ODE model. As before, capabilities of an EvoLudo module are advertised by implementing the appropriate interfaces. In order to advertise the ODE model, we need to implement the HasDE.ODE interface:

together with the accompanying import statement:

The interface HasDE.ODE requires implementing one additional method, which is easily accomplished in VS Code by hovering over the error location and selecting Quick fix…, followed by Add unimplemented methods:

In modules where the overall population density remains constant and hence the dynamics can be expressed in terms of frequencies, the method getDependent() returns the index of the dependent trait. In principle, this can be any one of the three traits $S$, $I$, or $R$. However, because the investigations will always focus on the infected or susceptible compartments, $I$ or $S$, let’s choose $R$ as the dependent trait.

In case the $R$ compartment is not relevant for the model, that is, if the $R$ trait is disabled, or no recruitment into the $R$ compartment happens, $p_{IR}=0$, the dependent trait should be $S$, keeping the focus on the infected trait $I$:

In order to implement the dynamical equations, Eq. (14.2), we need to define four additional fields that determine the transition rates between the different compartments, and assign them meaningful default values:

The exact values are not important but with those above, the infection becomes endemic and the disease persists in the population for most initial configurations.

What are the equilibrium frequencies of the three compartments? You can either solve Eq. (14.2) analytically by setting the left-hand side to zero and solve for the equilibria $S^{\ast },I^{\ast },R^{\ast }$ – or you can use EvoLudo to find the equilibria numerically. However, a few more steps are necessary before our module can accomplish this goal.

First and foremost we need to implement Eq. (14.2). This is done by defining a new inner class, say ODE, which extends the class org.evoludo.simulator.models.ODE.

More precisely, the class org.evoludo.simulator.models.ODE implements a simple numerical integrator based on Euler’s method. However, we can do significantly better by extending the RungeKutta class, which uses the more sophisticated fifth-order Cash-Karp Runge-Kutta method with variable step size (Press et al., 2007).

Plus the corresponding import statement for the RungeKutta class:

As before, we can take advantage of VS Code’s Quick fix… feature to create a stub for the constructor:

Since nothing else needs to be done in the constructor, the TODO comment has already been deleted in the code snippet above.

By default, the only view available in EvoLudo is the Console (see section 11.1.2: Views), which reports milestone events in the life cycle of a module (see section 12.10.1: Milestone listeners) as well as warnings and errors. Adding more detailed views is often as simple as implementing a particular interface.

With the successful implementation of the ODE model, adding a stochastic version of the $SIR$ model is rather straight forward. First we need to implement the HasDE.SDE interface to advertise to EvoLudo that the SIRE module also provides a model based on stochastic differential equations:

Just as for the ODE model, we now need to provide our own implementation of the SDE model. The part that is specific to the SIRE module are the rates of change that correspond to Eq. (14.2). Thus, let’s use the ODE model as a template and simply duplicate the inner class SIRE.ODE. Of course VS Code now complains about a duplicate nested type (SIRE.ODE). Simply replace all occurrences of ODE by SDE in the duplicated section to obtain:

Next, we need to return an instance of the inner class SIRE.SDE in the method createModel(Type type) by simply adding another if-condition:

Compile the module by reloading the EvoLudo page in your browser and verify that a new model is now available by clicking on Settings and Help. Now change to the SDE model by changing the EvoLudo settings to --model SDE and run the model. What are your observations?

Oops, quite clearly, the results are identical to the ODE model… What happened? Well, on closer inspection this is not surprising because we simply copied the ODE model!

Well, the implementations of the ODE and SDE models are almost identical. In fact, the only difference is the class that they extend but the key part of the code, namely the method getDerivatives(...) is identical in both cases. This is bad for code maintenance because for any changes to the ODE model, we have to remember to change the SDE model as well, and vice versa. This is not only tedious but also error-prone.

The easiest way to address this is to copy the method getDerivatives(...) to the outer class SIRE because then the method is accessible to both the ODE and SDE models. Also, while we are in the process of simplifying things, we can remove the double[] unused argument because it’s, well, unused, but add an argument accuracy because this field is unknown to the outer class. The modified method is:

With the streamlined ODE and SDE models in place, adding a PDE model is a piece of cake. First we need to implement the HasDE.PDE interface:

Next we need to provide our own implementation of the PDE model, SIRE.PDE. To start, let’s duplicate the code of the SDE model and rename all occurrences of SDE to PDE:

Finally, the method createModel(Type type) must return an instance of the inner class SIRE.PDE for the type Type.PDE by simply adding another if-condition:

That’s it! Done. However, what’s missing is a proper way to visualize the spatial dynamics and the distribution of the different types.

As always, when adding new features to an EvoLudo module, we start by implementing the appropriate interface. In this case, we want to add an individual based simulation model, IBS, and hence implement the HasIBS interface:

The default implementation of individual based simulation models in org.evoludo.simulator.models.IBS works fine for our purposes, and so we simply delegate the model creation to the parental class in the method createModel(Type type):

However, we do want to customize the interactions between individuals in the population. This effectively defines the disease progression in the $SIR$ model. The IBS can handle multiple species – a feature we will cover in greater detail in the next tutorial, chapter LABEL:chap:lv: LABEL:chap:lv. For each species IBS creates a population of individuals by calling the method createIBSPopulation(). This provides an opportunity for subclasses to return their custom implementation. In order to do just that, we create a subclass of IBSDPopulation, where the D refers to the fact that we are dealing with a discrete set of types, and return an instance of it in the method createIBSPopulation():

together with the import statement at the top:

Similar to the getDerivates(...) method in differential equation models, the core method in individual based simulation models is the method updatePlayerAt(int me, int[] refGroup, int rGroupSize). It gets called in every microscopic update step to update the trait of the focal individual, with the following arguments:

me:

the index of the focal individual,

refGroup:

the array of indices of the reference group members,

.

.

Attention: Never change the contents of the refGroup array! Depending on the parameter settings this may destroy the population geometry.

.

rGroupSize:

the size of the reference group.

.

.

Note: Usually refGroup.length will be different from rGroupSize. This allows to reuse the reference group array and reduce memory allocation.

.

return:

true if the focal individual updated its trait, false otherwise.

.

.

Note: A return value of true does not necessarily indicate an actual change of the focal individual’s type, just that the type of a model individual has been adopted.

.

The selection of the focal individual as well as of the reference group is controlled by various parameter settings, which includes the population update type, or the structure of the population. For a comprehensive list of parameters, see chapter 11: EvoLudo Reference. Most importantly, at this point those details do not need to concern us, EvoLudo takes care of it for us. Instead, we can focus on the actual disease dynamics. The focal individual is the only one that may update its trait in each microscopic update step. Thus, depending on the trait of the focal individual, different events can take place:

1. (i)

   A susceptible individual can become infected, provided the reference group contains at least one infected individual;

2. (ii)

   An infected individual can recover and become resistant, or it can return to the susceptible compartment;

3. (iii)

   A resistant individual may lose its immunity and return to the susceptible compartment.

The following code implements this logic:

plus the import statement:

Periodically the IBS model checks whether the simulation has converged, by calling the checkConvergence() method. Unlike in differential equation models, there is no general recipe to determine convergence in simulations. In our case, the simulation has converged if everyone is susceptible, i.e. the infection has cleared, or any mixture of susceptible and resistant individuals, if resistance is permanent, $p_{RS}=0$.

The addition of command line options for greater flexibility for exploring the dynamics of an EvoLudo module has already been covered in the previous tutorial, see section 13.4.1: Adding command-line options. For the SIRE module, we use the same option constructor but appropriately adjusted arguments. In order to set the infection rate, we define the following command-line option:

which requires several imports:

For individual based simulations, the parser ensures that the infection rate pSI is a probability. The other two options, for the rates of recovery and waning resistance, can be similarly implemented. Parsing the argument for --recover requires slightly more attention because of the optional second argument. The simplest approach is to parse the argument as a vector of doubles and then check the length of the vector.

An important extension of the basic $SIR$ model is to incorporate periodic variations in the transmission rates to model seasonal patterns, such as in influenza or Covid-19. This is easily accomplished by replacing the constant transmission rate, $p_{SI}$, with a periodically changing transmission function that depends on time $t$. This extension turns the homogeneous ODE in Eq. (14.2) into an inhomogeneous one. The simplest example is $p_{SI}(t)=p_{SI}+A\cdot \cos (2\pi \omega t)$, where $p_{SI}$ denotes the constant baseline infection rate, $A$ the amplitude and $\omega$ the frequency of the seasonal variation. For $A=0$ the constant transmission rate is recovered.

Thus, we require two additional parameters. Because they all relate to the transmission rate, let’s replace the double variable pSI by an array of three doubles with [pSI, A, ω]:

Hopf bifurcations and Turing patterns are fascinating dynamical phenomena that arise when fixed points lose their stability: on the one hand, when a stable focus loses its stability, a Hopf bifurcation can give rise to limit cycles and oscillatory behaviour; on the other hand, Turing patterns emerge when a stable, homogeneous equilibrium becomes unstable due to diffusion-driven instabilities, which can give rise to the spontaneous formation of spatial patterns. Turing instabilities are relevant for all kinds of pattern formation in nature, ranging from patchy vegetation in arid areas (Rietkerk et al., 2004, Scanlon et al., 2007) to the spots of the leopard (Gell-Mann, 1994). Both phenomena are of great interest in various fields, including epidemiology, ecology, and evolution.

Neither Hopf bifurcations nor Turing patterns are observed in our basic $SIR$ model. Depending on the parameters, the dynamics either converges to an endemic equilibrium, $Q$ – if it exists, it is globally stable – or the dynamics ends up in the disease-free equilibrium. However, it is relatively straightforward to modify the model to exhibit both Hopf bifurcations and Turing patterns by introducing non-linearities in the infection process. For example, by considering non-linear incidence rates, $p_{SI}{I}^p{S}^q$, which recovers our standard bi-linear setup for $p=q=1$ (Liu et al., 1986).

Because our population size is conserved, we have actually only two dynamical variables, say $S$ and $I$, while $R$ can be eliminated through $R=1-S-I$. Rewriting Eq. (14.2) accordingly, and including the non-linear incidence rates, we obtain:

Solving this dynamical system in general is a challenging task. However, at least for integer $p$ and $q$ we can still get some analytical traction, as we will see.

Note that the primary motivation to introduce these non-linearities was to better match empirical observations of disease outbreaks. They were first considered by Wilson and Worcester (1945b, a), who emphasized that they made this change “to investigate the consequences of various assumptions when the laws are not known.” This is important because it implies that the mechanistic processes that would give rise to such non-linearities are not well understood. As a consequence, what amounts to a small change to the ODE model can be very challenging to implement in individual based simulations. We have seen in section 2.8: From finite to infinite populations and section 5.5: From finite to infinite populations that the mass-action principle, or a Kramers-Moyal system size expansion, provide a deterministic way to derive an ODE from microscopic reaction equations. However, the reverse process, i.e. to derive microscopic rules from a given ODE, is not unique in general: different microscopic rules can give rise to the same macroscopic dynamics.

As an example, one possibility to achieve a non-linear incidence rate, based on a mechanistic description, is to assume that a susceptible individual needs to come into contact with multiple infected individuals simultaneously to contract the disease. For two infected individuals, the reaction kinetics can be written as a three-body interaction:

and results in an incidence rate of $p_{SI}S{I}^2$. Arguments along these lines can only justify integer exponents $p$ or $q$, while their biological, or epidemiological, interpretation and justification still needs to be supplied. Here we simply circumvent this problem by focussing solely on differential equations, and stick to the bi-linear incidence rates for individual based simulations.

In the context of the $SIR$ model, Hopf bifurcations can describe periodic outbreaks of a disease even in the absence of seasonal forcing (c.f. section 14.9: Extension I: Adding seasonal variation), while Turing patterns can spontaneously give rise to stable, but spatially heterogeneous, distributions of susceptible, infected, and recovered individuals, in cases where the ODE model predicts a stable, endemic equilibrium, $Q$.

In order to implement non-linear incidence rates, we introduce a new array incidence, which holds the two additional parameters, the exponents $p$ and $q$, and defaults of $1$ for each: