Practicum 02 - Model analysis, studying large-scale models in R

5 April 2022

Julien Arino

Department of Mathematics & Data Science Nexus
University of Manitoba*

Canadian Centre for Disease Modelling
Canadian COVID-19 Mathematical Modelling Task Force
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.

Outline

  • Steps of the analysis
  • The basic stuff (well-posedness, DFE)
  • Epidemic models
  • Endemic models
  • Numerical investigations of large-scale systems
  • Compound matrices

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

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

Numerical investigations of large-scale systems

Not very difficult

  • As for the mathematical analysis: if you do things carefully and think about things a bit, numerics are not hard. Well: not harder than numerics in low-D
  • Exploit vector structure

Set up parameters

pop = c(34.017, 1348.932, 1224.614, 173.593, 93.261) * 1e+06
countries = c("Canada", "China", "India", "Pakistan", "Philippines")
T = matrix(data = 
             c(0, 1268, 900, 489, 200, 
               1274, 0, 678, 859, 150, 
               985, 703, 0, 148, 58, 
               515, 893, 144, 0, 9, 
               209, 174, 90, 2, 0), 
           nrow = 5, ncol = 5, byrow = TRUE)

Work out movement matrix

p = list()
# Use the approximation explained in Arino & Portet (JMB 2015)
p$M = mat.or.vec(nr = dim(T)[1], nc = dim(T)[2])
for (from in 1:5) {
  for (to in 1:5) {
    p$M[to, from] = -log(1 - T[from, to]/pop[from])
  }
  p$M[from, from] = 0
}
p$M = p$M - diag(colSums(p$M))

p$P = dim(p$M)[1]
p$eta = rep(0.3, p$P)
p$epsilon = rep((1/1.5), p$P)
p$pi = rep(0.7, p$P)
p$gammaI = rep((1/5), p$P)
p$gammaA = rep((1/3), p$P)
# The desired values for R_0
R_0 = rep(1.5, p$P)

Write down indices of the different state variable types

Save index of state variable types in state variables vector (we have to use a vector and thus, for instance, the name "S" needs to be defined)

p$idx_S = 1:p$P
p$idx_L = (p$P+1):(2*p$P)
p$idx_I = (2*p$P+1):(3*p$P)
p$idx_A = (3*p$P+1):(4*p$P)
p$idx_R = (4*p$P+1):(5*p$P)

Set up IC and time

# Set initial conditions. For example, we start with 2
# infectious individuals in Canada.
L0 = mat.or.vec(p$P, 1)
I0 = mat.or.vec(p$P, 1)
A0 = mat.or.vec(p$P, 1)
R0 = mat.or.vec(p$P, 1)
I0[1] = 2
S0 = pop - (L0 + I0 + A0 + R0)
# Vector of initial conditions to be passed to ODE solver.
IC = c(S = S0, L = L0, I = I0, A = A0, R = R0)
# Time span of the simulation (5 years here)
tspan = seq(from = 0, to = 5 * 365.25, by = 0.1)

Set up to avoid blow up

Let us take for patches in isolation. Solve for 

for (i in 1:p$P) {
  p$beta[i] = 
    R_0[i] / S0[i] * 1/((1 - p$pi[i])/p$gammaI[i] + p$pi[i] * p$eta[i]/p$gammaA[i])
}

Define the vector field

SLIAR_metapop_rhs <- function(t, x, p) {
  with(as.list(p), {
    S = x[idx_S]
    L = x[idx_L]
    I = x[idx_I]
    A = x[idx_A]
    R = x[idx_R]
    N = S + L + I + A + R
    Phi = beta * S * (I + eta * A) / N
    dS = - Phi + MS %*% S
    dL = Phi - epsilon * L + p$ML %*% L
    dI = (1 - pi) * epsilon * L - gammaI * I + MI %*% I
    dA = pi * epsilon * L - gammaA * A + MA %*% A
    dR = gammaI * I + gammaA * A + MR %*% R
    dx = list(c(dS, dL, dI, dA, dR))
    return(dx)
  })
}

And now call the solver

# Call the ODE solver
sol <- ode(y = IC, 
           times = tspan, 
           func = SLIAR_metapop_rhs, 
           parms = p,
           method = "ode45")

One little trick (case with demography)

Suppose demographic EP is

Want to maintain for all to ignore convergence to demographic EP. Think in terms of :

So take

Then

and thus if , then and thus for all , i.e., for all

Word of warning about that trick, though..

has nonnegative (typically positive) diagonal entries and nonpositive off-diagonal entries

Easy to think of situations where the diagonal will be dominated by the off-diagonal, so could have negative entries

use this for numerics, not for the mathematical analysis

Compound matrices

The compound matrix method

  • An extension of Dulac's criterion to higher order systems
  • Useful to rule out the existence of periodic orbits
  • Was very popular for a while, but you must be aware of the main limitation:
    • Becomes hard to use when dimensionality

Refer for example to Fiedler for details

Let , , an -matrix and an integer,

Let and , and the lexicographically ordered sets of ordered -tuples of elements of and , respectively

The th compound matrix of is the matrix, with rows indexed by elements of and columns indexed by the elements of , s.t. element , , is the determinant

is the submatrix of consisting of rows in and columns in

Another interpretation of is as the th exterior product of matrix

Suppose now that is an -matrix. The matrix is a -matrix, with each entry a polynomial of with degree at most

Grouping monomials of same degree in

where matrices do not depend on

Matrix is the th additive compound matrix of and is denoted . It satisfies

The latter equality can be written

where is the right derivative

Suppose . Then, for

where is the entry in , is the entry in and is the number of entries of between and

When , we have

Suppose that . For all , let the th element of the lexicographic order of pairs of integers s.t.

Then the entry in the th row and th column of is

where is the entry in , is the entry in and is the number of elements of between and

Example

Let

Then

Let be two -matrices. Then

  • The number of nonzero offdiagonal entries of equals times the number of nonzero offdiagonal entries of
  • ,
  • (whence the additive suffix)
  • Let be a nonsingular -matrix. Then

Let be an -matrix. Denote the second additive compound matrix of

Let be a real -mix. For to have all its eigenvalues with negative real parts, it is necessary and sufficient that

  1. the second additive compound matrix has all its eigevalues with negative real parts

Role in stability

Consider the differential equation

A sufficient condition for a periodic orbit of be asymptotically orbitally stable with asymptotic phase is that the linear system

be asymptotically stable