Skip to contents

SquIPz

SquIPz.Rmd

SquIPz is a dynamical system describing the dynamics of the multiplicity of infection (MoI) in a cohort of humans as it ages. SquIPz was developed by extending SquIPz to include heterogeneous exposure as a Tweedie Process. In-Line documentation for SquIPz is accessed as:

help(SquIPz)

Implementation — The file SquIPz.R defines two functions:

dSquIPz 
solve_SquIPz
make_Z

This numerical implementation truncates the infinite system of differential equations. To verify the code and check the accuracy of the truncated system, we also derive and compute hybrid variables describing the first few moments of the infinite system.

Related

The Model

Variables

The independent variable is age, a.a.

The dependent variables are:

  • SS is the expected number uninfected and susceptible to infection

  • PP is the expected number uninfected and chemo-protected

  • IiI_i is the expected number infected with ii clones

We let HH denote the expected population density of the cohort as it ages: H=S+P+∑iIi.H = S+P+\sum_i I_i.

Similarly, we introduce the variable II:

I=∑iIiI = \sum_i I_i

Initial Conditions

In the model, we assume that everyone is born susceptible. HH is passed as a parameter, and the initial conditions are set to

  • S(0)=HS(0) = H

  • P(0)=0P(0) = 0

  • Ii(0)=0I_i(0) = 0

Mortality

The model assumes that all individuals die as they age at the rate μ\mu, such that

dHda=−μH\frac{\textstyle{dH}}{\textstyle{da}} = -\mu H

Exposure and Infection

We let hh denote the force of infection (FoI). This model, SquIPz, assumes that each incident infection increases the MoI by exactly one. The model SquIPz implements exposure as a Tweedie process.

The model assumes that each clone clears independently of the others at rate r,r, so for those infected with ii clones, the MoI goes down by one at the rate ri.ri.

Treatment and Chemprotection

These equations assume that individuals are treated and cured for several different reasons:

  • Incident infections cause disease and a fraction ρ\rho gets treated

  • Everyone in the population takes drugs at a background rate ξ\xi

  • In addition to other modes of treatment, infected individuals get treated at the higher rate σ\sigma

After being treated and cured, individuals enter the chemoprotected class P.P. Chemoprotection is lost at the rate η,\eta, and they enter the susceptible class S.S.

Differential Equations:

Changes in the susceptible population are:

dSda=ηP+rI1−(h+ξ+μ)S\frac{\textstyle{dS}}{\textstyle{da}} = \eta P + r I_1 - (h + \xi + \mu) S

Changes in the chemoprotected population are:

dPda=(hρ+ξ)(H−P)+σ∑iIi−(η+μ)P\frac{\textstyle{dP}}{\textstyle{da}} = \left(h \rho + \xi \right) \left(H-P \right) + \sigma \sum_i I_i - \left(\eta + \mu \right) P

Changes in the infected population are given by an infinite system of equations. The first one is a special case:

dI1da=h(1−ρ)S+2rI2−(r+h+ξ+σ+μ)I1\frac{\textstyle{dI_1}}{\textstyle{da}} = h(1-\rho)S + 2 r I_2 - (r + h + \xi + \sigma + \mu) I_1

for i>1i>1:

dIida=h(1−ρ)Ii−1+(i+1)rIi+1−(ir+h+ξ+σ+μ)I1\frac{\textstyle{dI_i}}{\textstyle{da}} = h(1-\rho)I_{i-1} + (i+1) r I_{i+1} - (i r + h + \xi + \sigma + \mu) I_1

Implementation and Verification

In this implementation, the infinite system of equations is truncated at N.N. For the last differential equation IN,I_N, no new infections occur (h=0h=0) and there is no IN+1I_{N+1} so nothing is added from above:

dINda=h(1−ρ)IN−1−(Nr+ξ+σ+μ)IN\frac{\textstyle{dI_N}}{\textstyle{da}} = h(1-\rho)I_{N-1} - (N r + \xi + \sigma + \mu) I_N

In the implementation, if no value for NN is passed, it is set to

N = round(max(10*h/r,20))

Hybrid Variables:

This implementation also computes several other variables: Let I=∑iIiI = \sum_i I_i

The dynamics of the hybrid variable II are:

dIda=h(1−ρ)S−rI1−(ρh+ξ+σ+μ)I\frac{\textstyle{dI}}{\textstyle{da}} = h(1-\rho)S - r I_1 - (\rho h+\xi+\sigma+\mu)I

Let mnm_n denote the nthn^{th} moment of the distribution of the MoI:

mn=1H∑iinIim_n = \frac{\textstyle{1}}{\textstyle{H}}\sum_i i^n I_i

With some algebra (see below), we can derive an equation for the dynamics of the hybrid variable m1m_1 are: dm1da=h(1−ρ)(1−PH)−(r+hρ+σ+ξ)m1\frac{\textstyle{dm_1}}{\textstyle{da}} = h(1-\rho)\left(1-\frac{\textstyle{P}}{\textstyle{H}}\right) - \left(r + h \rho +\sigma + \xi\right) m_1

Similarly, we can compute the dynamics of m2:m_2:

dm2da=h(1−ρ)(1−PH+2m1)−(2r+hρ+σ+ξ)m2+rm1\frac{\textstyle{dm_2}}{\textstyle{da}} = h(1-\rho)\left( 1- \frac{\textstyle{P}}{\textstyle{H}} + 2 m_1\right) - \left(2r + h \rho + \sigma + \xi \right) m_2 + r m_1

The dynamics of the hybrid variables, I,I,m1,m_1, and m2m_2 are verified by solving the equations and using the formulas to compute their values.

Derivations

Change of Variables

I=∑iIim=∑iiIi/Hm2=∑ii2Ii/H. \begin{array}{ccc} I &= & \sum_i I_i \\[6pt] m &= & \sum_i i I_i /H\\[6pt] m_2 &= & \sum_i i^2 I_i /H\\ \end{array}.

Infected

dIda=h(1−ρ)S∑izi−∑i(h+σ+ξ+μ)Ii−rI1−∑i−riIi+(i+1)rIi+1+h(1−ρ)∑i∑j<iIi−jzj\frac{dI}{da} =\begin{array}{l} h(1-\rho) S \sum_i z_i \\ - \sum_i \left(h +\sigma + \xi + \mu\right) I_i \\ - r I_1 - \sum_i -ri I_i + (i+1) r I_{i+1} \\ + h(1-\rho)\sum_i \sum_{j<i} I_{i-j} z_j\\ \end{array} We can use some identities to simplify a lot:

dIda=h(1−ρ)S−(h+σ+ξ+μ)I−r1I1−h(1−ρ)∑i∑j<iIi−jzj\frac{dI}{da} =\begin{array}{l} h(1-\rho) S - \left(h + \sigma + \xi + \mu\right) I - r_1 I_1\\ - h(1-\rho)\sum_i \sum_{j<i} I_{i-j} z_j\\ \end{array} but ∑i∑j<iIi−jzj=∑i∑jzjIi=I\sum_i \sum_{j<i} I_{i-j} z_j = \sum_i \sum_j z_j I_i = I so dIda=h(1−ρ)S−(hρ+σ+ξ+μ)I−r1I1\frac{dI}{da} = h(1-\rho) S - \left(h \rho + \sigma + \xi + \mu\right) I - r_1 I_1

Mean MoI

dmda=h(1−ρ)S∑iizi/H−∑i(h+σ+ξ+μ)iIi/H∑ii((i+1)rIi+1−riIi)/H+h(1−ρ)∑i∑j<iiIi−jzj\frac{dm}{da} =\begin{array}{l} h(1-\rho) S \sum_i i z_i/H\\ - \sum_i \left(h +\sigma + \xi + \mu\right) i I_i /H \\ \sum_i i \left((i+1)r I_{i+1} - r_i I_i \right) /H \\ + h(1-\rho)\sum_i \sum_{j<i} i I_{i-j} z_j\\ \end{array}

rr terms

For the terms involving rr, we want to simplify ∑ii((i+1)rIi+1−riIi) \sum_i i \left((i+1) r I_{i+1} - r i I_i \right) In long form, this is just 1(−rI1+2rI2)2(−2rI2+3rI3)3(−3rI3+4rI4) \begin{array}{l} 1\left(-r I_1 + 2r I_2 \right) \\ 2 \left(-2r I_2 + 3r I_3\right) \\ 3 \left(-3r I_3 + 4r I_4\right) \\ \end{array} but this simplifies to −rI1−2rI2−3rI3−…=−rm-r I_1 - 2r I_2 - 3r I_3 - \ldots = -rm

dmda=h(1−ρ)Sẑ/H−∑i(h+σ+ξ+μ)m−rI1+∑i((i+1)r(i+1)Ii+1−riIi)/H+h(1−ρ)∑i∑j<iiIi−jzj\frac{dm}{da} =\begin{array}{l} h(1-\rho) S \hat z/H\\ - \sum_i \left(h +\sigma + \xi + \mu\right) m \\ - r I_1 + \sum_i \left((i+1) r (i+1) I_{i+1} - r_i I_i \right) /H \\ + h(1-\rho)\sum_i \sum_{j<i} i I_{i-j} z_j\\ \end{array}

hh terms

h(1−ρ)∑i∑j<ii2Ii−jzjh(1-\rho)\sum_i \sum_{j<i} i^2 I_{i-j} z_j In long form, summing from IiI_i we get I1−hI1+h(1−ρ)I1(2z1+3z2+4z3+…)I2−22hI2+h(1−ρ)I2(32z1+42z2+52z3+…)⋮⋮Ii−i2hIi+h(1−ρ)Ii((i+1)2z1+(i+2)2z2+…)\begin{array}{c|l} I_1 & - h I_1 + h (1-\rho) I_1 \left(2 z_1 + 3 z_2 + 4 z_3 + \ldots\right) \\ I_2 & - 2^2 h I_2 + h (1-\rho) I_2 \left(3^2 z_1 + 4^2 z_2 + 5^2 z_3 + \ldots\right) \\ \vdots & \vdots \\ I_i & - i^2 h I_i + h (1-\rho) I_i \left((i+1)^2 z_1 + (i+2)^2 z_2 + \ldots \right) \end{array} but since
∑j(i+j)zj=∑jizj+∑jjzj=i+ẑ \sum_j (i+j) z_j = \sum_j i z_j + \sum_j j z_j = i + \hat z

We can rewrite all this above as: I1−hI1+h(1−ρ)I1(1+ẑ)I2−2hI2+h(1−ρ)I2(2+ẑ)⋮⋮Ii−ihIi+h(1−ρ)Ii(i+ẑ)\begin{array}{c|l} I_1 & - h I_1 + h (1-\rho) I_1 \left(1 + \hat z\right) \\ I_2 & - 2 h I_2 + h (1-\rho) I_2 \left(2 + \hat z\right) \\ \vdots & \vdots \\ I_i & - i h I_i + h (1-\rho) I_i \left(i + \hat z \right) \end{array} now dividing by H,H, we get dmda=−hm+h(1−ρ)m+h(1−ρ)ẑS+IH+terms not involving h \frac{dm}{da} = - h m + h (1-\rho)m + h(1-\rho)\hat z \frac{S+I}{H} + \mbox{terms not involving } h

MoI, Second Moment

dmda=h(1−ρ)S∑ii2zi/H−∑i(h+σ+ξ+μ)i2Ii/H∑ii2((i+1)rIi+1−riIi)/H+h(1−ρ)∑i∑j<ii2Ii−jzj\frac{dm}{da} =\begin{array}{l} h(1-\rho) S \sum_i i^2 z_i/H\\ - \sum_i \left(h +\sigma + \xi + \mu\right) i^2 I_i /H \\ \sum_i i^2 \left((i+1)r I_{i+1} - r_i I_i \right) /H \\ + h(1-\rho)\sum_i \sum_{j<i} i^2 I_{i-j} z_j\\ \end{array}

hh terms

h(1−ρ)∑i∑j<iiIi−jzjh(1-\rho)\sum_i \sum_{j<i} i I_{i-j} z_j In long form, summing from IiI_i we get I1−hI1+h(1−ρ)I1(22z1+32z2+42z3+…)I2−22hI2+h(1−ρ)I2(32z1+42z2+52z3+…)⋮⋮Ii−i2hIi+h(1−ρ)Ii((i+1)2z1+(i+2)2z2+…)\begin{array}{c|l} I_1 & - h I_1 + h (1-\rho) I_1 \left(2^2 z_1 + 3^2 z_2 + 4^2 z_3 + \ldots\right) \\ I_2 & - 2^2 h I_2 + h (1-\rho) I_2 \left(3^2 z_1 + 4^2 z_2 + 5^2 z_3 + \ldots\right) \\ \vdots & \vdots \\ I_i & - i^2 h I_i + h (1-\rho) I_i \left((i+1)^2z_1 + (i+2)^2 z_2 + \ldots \right) \end{array} but since
∑j(i+j)2zj=∑ji2zj+∑j2ijzj+∑jj2zj=i2+2iẑ+ẑ2 \sum_j (i+j)^2 z_j = \sum_j i^2 z_j + \sum_j 2 i j z_j + \sum_j j^2 z_j = i^2 + 2 i \hat z + \hat z_2

We can rewrite all this above as: I1−hI1+h(1−ρ)I1(1+1⋅2ẑ+ẑ2)I2−22hI2+h(1−ρ)I2(22+2⋅2ẑ+ẑ2)⋮⋮Ii−i2hIi+h(1−ρ)Ii(i2+i⋅2ẑ+ẑ2)\begin{array}{c|l} I_1 & - h I_1 + h (1-\rho) I_1 \left(1 + 1 \cdot 2 \hat z + \hat z_2\right) \\ I_2 & - 2^2 h I_2 + h (1-\rho) I_2 \left(2^2 + 2 \cdot 2 \hat z + \hat z_2 \right) \\ \vdots & \vdots \\ I_i & - i^2 h I_i + h (1-\rho) I_i \left(i^2 + i \cdot 2 \hat z + \hat z_2 \right) \end{array} now dividing by H,H, we get dm2da=−hm2+h(1−ρ)(m2+2ẑm+ẑ2IH)+h(1−ρ)ẑ2SH+terms not involving h \frac{dm_2}{da} = - h m_2 + h (1-\rho)\left(m_2 + 2 \hat z m + \hat z_2 \frac I H \right) + h (1-\rho) \hat z_2 \frac{S}{H} + \mbox{terms not involving } h

The old correct formula was: dm2da=h(1−ρ)ẑ(S+IH+2m)−(2r+hρ+σ+ξ)m2+rm\frac{dm_2}{da} = h(1-\rho)\hat z \left(\frac{S + I}H + 2 m \right) - (2r + h\rho + \sigma + \xi)m_2 + r m

Numerical Verification 

F_z = make_NBtrunc(1, 2)
plot(F_z(1:10), type = "h")

F_a = make_F_a(1)
q_out <- solve_SquIPz(8/365, F_a, F_z, sigma=2/365, xi=0, Amax=3*365, N=400)

For verification, the variable m1m_1 is computed two ways, and we
we note that they match. To see it, we plotted both, one in solid black and the other dashed yellow:

with(q_out, plot(age, m_1, type = "l", ylab = expression(m[1]), main = "First Moment of the MoI"))
with(q_out, lines(age, m1, lty=2, col = "yellow"))

Since it’s hard to see differences, we can simply plot the differences. Here, we’ve set the limits to be 10−1410^{-14} so we can visualize the differences:

ylm = 1e-13
with(q_out, plot(age, m_1-m1, type = "l", ylab = "Errors", main = "First Moment: Numerical Errors", ylim = c(-ylm, ylm)))

with(q_out, plot(age/365, m_2, type = "l", main = "Second Moment of the MoI",
                 ylab = expression(m[2]), 
                 xlab = "Age", ylim = range(m_2, m2)))
with(q_out, lines(age/365, m2, lty=2, col = "yellow"))

ylm = 1e-10/5
with(q_out, plot(age, m_2-m2, type = "l", 
                 ylab = "Errors", 
                 main = "Second Moment: Numerical Errors", 
                 ylim = c(-ylm, ylm)))

We note that if we had truncated the system of equations at N=4,N=4, the moments diverge:

q1_out <- solve_SquIPz(8/365, F_a, F_z, sigma=2/365, xi=0, Amax=3*365, N=4)
with(q1_out, plot(age, m_1, type = "l", 
                  ylab = expression(m[1]), 
                  main = "First Moment of the MoI, N=4"))
with(q1_out, lines(age, m1, lty=2, col = "yellow"))

# Example

clrs = viridisLite::turbo(7)
set.seed(234)
Sa = makepar_F_type2()
Sp = makepar_F_sin()
Tp = makepar_F_spline(seq(0, 3650, length.out=11), 1+runif(11, -1, 1), X=2)
Kp = makepar_F_sharkbite(D=730, L=365)
F_t <- make_ts_function(scale = 0.05, season_par=Sp, trend_par=Tp, shock_par=Kp)
F_a <- make_F_a(avg = 3/365, age_par=Sa, season_par=Sp, trend_par=Tp, shock_par=Kp)

tt <- seq(0, 3650, by=5)
aa <-seq(0, 365*5, by =5) 
plot(tt, 0.05*F_t(tt), type = "l")

plot(tt, 0.1*F_t(tt), type = "l", lwd=2, ylim = c(0,0.04))
lines(aa, F_a(aa), type = "l", col = clrs[2])
lines(aa+365, F_a(aa, 365), type = "l", col = clrs[3])
lines(aa+730, F_a(aa, 730), type = "l", col = clrs[4])
lines(aa+1095, F_a(aa, 1095), type = "l", col = clrs[5])
lines(aa+1460, F_a(aa, 1460), type = "l", col = clrs[6])

plot(aa, F_a(aa), type = "n", lwd=2, ylim = c(0,0.04))
lines(aa, F_a(aa), type = "l", col = clrs[2])
lines(aa, F_a(aa, 365), type = "l", col = clrs[3])
lines(aa, F_a(aa, 730), type = "l", col = clrs[4])
lines(aa, F_a(aa, 1095), type = "l", col = clrs[5])
lines(aa, F_a(aa, 1460), type = "l", col = clrs[6])

plot(aa, cumsum(F_a(aa)), ylab = "Age", xlab = "Cumulative Exposure", type = "n", ylim = c(0,3))
lines(aa, cumsum(F_a(aa)), type = "l", col = clrs[2])
lines(aa, cumsum(F_a(aa, 365)), type = "l", col = clrs[3])
lines(aa, cumsum(F_a(aa, 730)), type = "l", col = clrs[4])
lines(aa, cumsum(F_a(aa, 1095)), type = "l", col = clrs[5])
lines(aa, cumsum(F_a(aa, 1460)), type = "l", col = clrs[6])

q3_out1 <- solve_SquIPz(5/365, F_a, F_z, sigma=2/365, xi=0, Amax=5*365)
q3_out2 <- solve_SquIPz(5/365, bday=365, F_a, F_z, sigma=2/365, xi=0, Amax=5*365)
q3_out3 <- solve_SquIPz(5/365, bday=730, F_a, F_z, sigma=2/365, xi=0, Amax=5*365)
q3_out4 <- solve_SquIPz(5/365, bday=365*3, F_a, F_z, sigma=2/365, xi=0, Amax=5*365)
q3_out5 <- solve_SquIPz(5/365, bday=365*4, F_a, F_z, sigma=2/365, xi=0, Amax=5*365)
with(q3_out1, plot(age, m_1/x, col = clrs[2], ylab = "mean MoI", type ="l", ylim = c(1, 2)))
with(q3_out2, lines(age, m_1/x, col = clrs[3]))
with(q3_out3, lines(age, m_1/x, col = clrs[4]))
with(q3_out4, lines(age, m_1/x, col = clrs[5]))
with(q3_out5, lines(age, m_1/x, col = clrs[6]))