Basic steps in an analysis of an epidemiological model

26 August 2023

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

Canadian Centre for Disease Modelling
PHAC-EMNID / REMMI-ASPC
NSERC-PHAC EID Modelling Consortium (CANMOD, MfPH, OMNI/RÉUNIS)

*The University of Manitoba campuses are located on original lands of Anishinaabeg, Cree, Oji-Cree, Dakota and Dene peoples, and on the homeland of the Métis Nation. We respect the Treaties that were made on these territories, we acknowledge the harms and mistakes of the past, and we dedicate ourselves to move forward in partnership with Indigenous communities in a spirit of reconciliation and collaboration.

Outline

  • Steps of the analysis
  • The basic stuff (well-posedness, DFE)
  • Epidemic models
  • Endemic models

Steps of the analysis

Step 0 - Well-posedness

  • Do solutions exist?
  • Are they unique?
  • Are they bounded?
  • Invariance of the nonnegative cone under the flow..?

For "classic" models, all of these properties are more or less a given, so good to bear in mind, worth mentioning in a paper, but not necessarily worth showing unless this is a MSc or PhD manuscript

When you start considering nonstandard models, or PDE/DDE, then often required

Step 1 - Epidemic model or endemic model ?

  • Often a source of confusion: analysis of epidemic models differs from analysis of endemic models!!
  • Important to determine what you are dealing with
  • Easy first test (can be wrong): is there demography?
    • Demography can lead to constant population, but if there is "flow" through the system (with, e.g., births = deaths), then there is demography
  • Other (more complex) test: what is the nature of the DFE?

Step 1 and a half - Computing the DFE

  • If you are not yet sure whether you have an epidemic or endemic model, you need to compute the DFE (you will need it/them anyway)
  • Usually: set all infected variables to 0 (I, L and I, etc.)
    • If you find a single or denumerable number of equilibria for the remaining variables, this is an endemic model
    • If you get something of the form "any value of works", this is an epidemic model

Step 2 - Epidemic case

  • Compute
  • Usually: do not consider LAS properties of DFE, they are given
  • Compute a final size (if feasible)

Step 2 - Endemic case

  • Compute and deduce LAS properties of DFE
  • (Optional) Determine direction of bifurcation at
  • (Sometimes impossible) Determine GAS properties of DFE or EEP

Why considering LAS properties of epidemic model is wrong

Consider the IVP

and denote its solution at time through the initial condition

is an equilibrium point if

is locally asymptotically stable (LAS) if open in the domain of s.t. for all , for all and furthermore,

If there is a continuum of equilibria, then , where is some curve in the domain of , s.t. for all . We say is not isolated. But then any open neighbourhood of contains elements of and taking , , implies that . is locally stable but not locally asymptotically stable!

The basic stuff (well-posedness, DFE)

Existence and uniqueness

  • Is your vector field ?
    • If so, you are done
    • If not, might be worth checking. Some of the models in particular have issues if the total population is variable and under circumstances
  • Probably not worth more than "solutions exist and are unique" in most instances...

Invariance of the nonnegative cone under the flow

  • Study of this can be warranted
  • What can be important is invariance of some subsets of the nonnegative cone under the flow of the system.. this can really help in some cases

Example: SIS system

First, remark that - is , giving existence and uniqueness of solutions

Invariance of under the flow of -

If , then becomes

cannot ever become zero: if , then for all . If , then for small and by the preceding argument, this is also true for all

For , remark that if , then is positively invariant: if , then for all

In practice, values of for any solution in are "carried" by one of the following 3 solutions:

  1. : increases to
  2. : remains equal to
  3. : decreases to

As a consequence, no solution with can enter . Suppose and s.t. ; denote the value of when becomes zero

Existence of contradicts uniqueness of solutions, since at , there are then two solutions: that initiated in and that initiated with

Boundedness

positive quadrant (positively) invariant under flow of -

We could detail more precisely (positive IC ..) but this suffices here

From the invariance dans the boundedness of the total population , we deduce that solutions to - are bounded

Where things can become complicated...

  • If , e.g., , what happens to the incidence?
  • If , e.g., , solutions are unbounded

Computing the DFE

  • Set all infected variables to zero, see what happens...
  • Personnally: I prefer to set some infected variables to zero and see if I recuperate the DFE that way

Epidemic models

  • Computation of
  • Final size relation
  • Examples

Arino, Brauer, PvdD, Watmough & Wu. A final size relation for epidemic models. Mathematical Biosciences and Engineering 4(2):159-175 (2007)

Computation of

A method for computing in epidemic models

  • This method is not universal! It works in a relatively large class of models, but not everywhere. If it doesn't work, the next generation matrix method (see later) does work, but should be considered only for obtaining the reproduction number, not to deduce LAS (cf. my remark earler)
  • Here, I change the notation in the paper, for convenience

Standard form of the system

Suppose system can be written in the form

where , and are susceptible, infected and removed compartments, respectively

IC are with at least one of the components of positive

  • continuous function encoding recruitment and death of uninfected individuals
  • diagonal with diagonal entries the relative susceptibilities of susceptible compartments, with convention that
  • Scalar valued function represents infectivity, with, e.g., for mass action
  • row vector of relative horizontal transmissions

  • has entry the fraction of individuals in susceptible compartment that enter infected compartment upon infection
  • diagonal with diagonal entries the relative susceptibilities of susceptible compartments, with convention that
  • Scalar valued function represents infectivity, with, e.g., for mass action
  • row vector of relative horizontal transmissions
  • describes transitions between infected states and removals from these states due to recovery or death

  • continuous function encoding flows into and out of removed compartments because of immunisation or similar processes
  • has entry the rate at which individuals in the infected compartment move into the removed compartment

Suppose is a locally stable disease-free equilibrium (DFE) of the system without disease, i.e., an EP of

Let

  • If , the DFE is a locally asymptotically stable EP of -
  • If , the DFE of - is unstable

If no demopgraphy (epidemic model), then just , of course

Final size relations

Final size relations

Assume no demography, then system should be writeable as

For continuous, define

Define the row vector

then

Suppose incidence is mass action, i.e., and

Then for , express as a function of using

then substitute into

which is a final size relation for the general system when

If incidence is mass action and (only one susceptible compartment), reduces to the KMK form

In the case of more general incidence functions, the final size relations are inequalities of the form, for ,

where is the initial total population

Examples

The SLIAR model

center

Here, , and , so , and

Incidence is mass action so and thus

For final size, since , we can use :

Suppose , then

If , then

A model with vaccination

Fraction of are vaccinated before the epidemic; vaccination reduces probability and duration of infection, infectiousness and reduces mortality

with and

Here, , ,

So

and the final size relation is

Where things go awry, final size-wise!

  • Summer 2021 work of Aaron Shalev (U of M)
  • Consider the very simple 2-patch metapopulation (see Lecture 05)

  • Unidirectional movement from patch 1 to patch 2
  • cannot be scalar-valued here!

Endemic models

  • Computation of using the next generation matrix method
  • Examples

Computation of using the next generation matrix method

Next generation matrix/operator

Diekmann and Heesterbeek, characterised in ODE case by PvdD & Watmough (2002)

Consider only compartments with infected individuals and write

  • flows into infected compartments because of new infections
  • other flows (with sign)

Compute the (Frechet) derivatives and with respect to the infected variables and evaluate at the DFE

Then

where is the spectral radius

Preliminary setup

, , with the first compartments the infected ones

the set of all disease free states:

Distinguish new infections from all other changes in population

  • rate of appearance of new infections in compartment
  • rate of transfer of individuals into compartment by all other means
  • rate of transfer of individuals out of compartment

Assume each function continuously differentiable at least twice in each variable

where

Some assumptions

  • (A1) If , then for

Since each function represents a directed transfer of individuals, all are non-negative

  • (A2) If then . In particular, if , then for

If a compartment is empty, there can be no transfer of individuals out of the compartment by death, infection, nor any other means

  • (A3) if

The incidence of infection for uninfected compartments is zero

  • (A4) If then and for

Assume that if the population is free of disease then the population will remain free of disease; i.e., there is no (density independent) immigration of infectives

One last assumption for the road

Let be a DFE of the system, i.e., a (locally asymptotically) stable equilibrium solution of the disease free model, i.e., the system restricted to . We need not assume that the model has a unique DFE

Let be the Jacobian matrix . Some derivatives are one sided, since is on the domain boundary

  • (A5) If is set to zero, then all eigenvalues of have negative real parts

Note: if the method ever fails to work, it is usually with (A5) that lies the problem

Stability of the DFE as function of

Suppose the DFE exists. Let then

with matrices and obtained as indicated. Assume conditions (A1) through (A5) hold. Then

  • if , then the DFE is LAS
  • if , the DFE is unstable

Important to stress local nature of stability that is deduced from this result. We will see later that even when , there can be several positive equilibria

Direction of the bifurcation at

bifurcation parameter s.t. for and for and DFE for all values of and consider the system

Write

as block matrix

Write , , the entry of and let and be left and right eigenvectors of s.t.

Let

Consider model with satisfying conditions (A1)–(A5) and as described above

Assume that the zero eigenvalue of is simple

Define and by and ; assume that . Then s.t.

  • if , then there are LAS endemic equilibria near for
  • if , then there are unstable endemic equilibria near for

Examples

  • The SLIRS model
  • A tuberculosis model incorporating treatment
  • An "issue" with the next generation method for
  • Bistable states

The SLIRS model

Example of the SLIRS model

Variations of the infected variables described by

Thus

Then compute the Jacobian matrices of vectors and

where

We have

Also, when constant, , then

and thus,

A tuberculosis model incorporating treatment

  • While undergoing treatment, individuals can be infected

C & Feng. To treat or not to treat: the case of tuberculosis. Journal of Mathematical Biology 35: 629-656 (1997).

with

DFE is with solution to (so there must hold that )

Assume without loss of generality that

So

and

Effect of another choice of (1)

So

and

Effect of another choice of (2)

So

and

has a simple zero eigenvalue when . All second derivatives of are zero at the DFE except

So

where the eigenvectors and can be chosen so that and

Typically, so bifurcation is forward

A slight change to the model

Feng, C & Capurro. A Model for Tuberculosis with Exogenous Reinfection. Theoretical Population Biology 57(3): 235-247 (2000)

With this change

where it can be shown that

Therefore, there are situations when , i.e., there can be backward bifurcations

An "issue" with the next generation method for

(Malaria model with transgenic vectors)

JA, Bowman, Gumel & Portet. Effect of pathogen-resistant vectors on the transmission dynamics of a vector-borne disease. Journal of Biological Dynamics 1:320-346 (2007)

Foreword

  • This is not an issue with the method itself, but an illustration of the reason why it is important to check that (A1)-(A5) are satisfied when using the next generation matrix method of PvdD & Watmough (2002)

Hypotheses about mosquitoes

  • H1 Resistant vectors are completely immune to the pathogen
  • H2-a  A proportion of the offspring resulting from the interbreeding of wild and resistant vectors are resistant to the pathogen
  • H2-b  A proportion $ of the offspring resulting from inbreeding within the resistant type are resistant to the pathogen
  • H3 In the absence of disease, wild vectors are better fitted for competition than resistant vectors
  • H4 Wild vectors are ill-fitted for competition when they carry the disease

The model

  • if or
  • if or

Coexistence of vectors

  • Can wild type and transgenic vectors coexist?

Fitness (from undetailed assumptions, ):

Find 2 boundary EP and and one coexistence EP . If there is no loss of resistance in offspring of two resistant parents (), then there is another boundary EP

Where the problem arises

There are up to 4 EPs for vectors and these are independent from the host population in the case when disease is absent

the whole coupled system with vectors and hosts has up to 4 DFE

We can use the method of PvdD & Watmough (2002) at each of these DFE to get the local stability properties of these DFE

At , we find , which makes no sense. What's wrong?

Problem is with (A5): compute Jacobian of system and evaluate at , we get e-values and , so is always unstable (A5) cannot be satisfied at and the LAS condition provided by PvdD & Watmough (2002) is not usable

Bistable states

Undesired effect of vaccination

Arino, McCluskey & PvdD. Global results for an epidemic model with vaccination that exhibits backward bifurcation. SIAM J Applied Math (2003)

Another SIRS model with vaccination

center

SIRSV epidemic model

  • proportion of newborns vaccinated
  • vaccination rate of susceptibles
  • rate of vaccine efficacy loss
  • vaccine efficacy

Since the total population is constant, the system in proportions takes the form

where , , , are the proportions of individuals who are in the susceptible, infectious, recovered and vaccinated, respectively

Equilibrium and bifurcations

The system always has the DFE

We now consider endemic equilibria with

When (vaccine 100% efficacious), there is at most one endemic equilibrium. From now on, assume (realistic) that , i.e., vaccine is not 100% efficacious

Existence of endemic equilbria

The existence of endemic equilibria is determined by the number of positive roots of the polynomial

where

Case of a forward bifurcation

Case of a backward bifurcation

Bistable region

  • Concavity of the curve is fixed (since ), so a necessary condition for existence of two endemic equilibria is:
    • and
    • The roots of must be real

region of bistability is , and

Bifurcation in the plane

center

EEP

If there are such solutions to (potentielly a double root), EEP are

Using the next generation method, the reproduction number (with vaccination) is

where

and as a consequence

Stability - DFE

  • Using a theorem of PvdD & Watmough (2002), the DFE is
    • locally asymptotically stable for
    • unstable for
  • Furthermore, when , using as a Lyapunov function, it is easily shown the the DFE is globally asymptotically stable

Local stability - EEP

Linearising at the EEP

  • at the smaller , Jacobian matrix has negative trace and positive determinant one of the eigenvalues is positive and the lower bifurcating branch is unstable
  • On the upper branch, conclude from linearisation that there is either one or three eigenvalues with nonpositive real part stability is undetermined. From numerical investigations, the upper branch seems locally stable

Spectral abscissa at the EP

Spectral abscissa (maximum of real parts of eigenvalues) of the linearisation at the DFE and the 2 EEP, when varies

center

Global behaviour

Suppose that in the system, parameters satisfy

Then all positive semi-trajectories in , where

limit to a unique equilibrium point