Preliminary example #1

Constant increase of Infected compartment: \[\frac{\mathrm dI(t)}{\mathrm dt} = a,\quad I(0)=I_0,\] which is solved by moving \(\mathrm{d}t\) to the right hand side and re-writing the equation in the form: \[\mathrm dI = a\mathrm dt \quad \Rightarrow \quad \int \mathrm dI = \int a\>\mathrm dt.\] This gives: \(I(t) = a\,t+C\), where \(C\) is a constant of integration. It is defined by using initial condition: \(I(0)=I_0\), so that: \(I_0=a\times0+C\) and \(C=I_0\).

Finally we write: \[I(t) = I_0+a\,t.\]

library(dplyr)
library(magrittr)
library(ggplot2)

# Initial number of susceptibles: 100 people
I0 = 5
# Constant rate
a = 3
# We set a time to a sequence of 100 timepoint 
# between 0 and 5 weeks
df = data.frame(t = seq(0,5,length.out=100)) %>%
# then we impute the solution we obtained before 
  mutate(I = I0 + a*t)

ggplot(data=df,aes(x=t)) +
  geom_path(aes(y=I,color="Infected"),size=1.5) +
  coord_cartesian(ylim=c(0,27)) +
  theme_bw(base_size=14,base_family='Times') +
  labs(x="Time (weeks)",y="Number of infected people") + 
  guides(color=guide_legend(title=NULL)) +
  theme(legend.background = element_blank())

Preliminary example 2

A small modification of the previous example: \[\frac{\mathrm dI(t)}{\mathrm dt} = a\cdot t,\quad I(0)=I_0,\] that is solved in analogous way by moving \(\mathrm dt\) to the right hand side: \[\int \mathrm dI(t) = \int a\cdot t{\mathrm dt}\quad \Rightarrow \quad I(t) = \frac{at^2}2+C\]. The constant of integration is defined by the initial condition: \(I(0)=I_0=\frac{a\times0^2}2+C\), so that: \(C=I_0\) and: \[I(t)=I_0+\frac{at^2}2\].

# Here we modify the dataframe from a previous example just a bit
# Note that the power of two is written as "**" in R
df %>% 
  mutate(I = I0 + a*t**2/2) %>%
  ggplot(aes(x=t)) +
    geom_path(aes(y=I,color="Infected"),size=1.5) +
    coord_cartesian(ylim=c(0,52)) +
    theme_bw(base_size=14,base_family='Times') +
    labs(x="Time (weeks)",y="Number of infected people") + 
    guides(color=guide_legend(title=NULL)) +
    theme(legend.background = element_blank())

Example 1: Catalytic model

Dynamics: \[\frac{\mathrm dS(t)}{\mathrm dt} = -\lambda S(t),\quad \frac{\mathrm dI(t)}{\mathrm dt} = \lambda S(t),\] with initial conditions: \(S(0)=S_0\) and \(I(0)=0\).

First, we solve the equation for \(S(t)\): \[ \frac{\mathrm dS}S = -\lambda\mathrm dt\quad\Rightarrow\quad \mathrm d\ln S = -\lambda\mathrm dt\,. \] ⊕Because: \(\mathrm d\ln S=\frac{\mathrm dS}S\)

Integrating both sides gives: \[ \int\mathrm d\ln S = -\lambda\int\mathrm dt\quad\Rightarrow\quad \ln S = -\lambda t+C_1\,, \] where \(C_1\) is the constant of integration.

Now we take the exponent of both sides: \[ \exp(\ln S) = \exp(-\lambda t+C_1)=\exp(-\lambda t)\cdot\exp(C_1)\,, \] or \[ S(t) = e^{C_1}\cdot e^{-\lambda t}\,. \] Here we do not know the constant \(C_1\), so we define it using initial condition. Precisely, we substitute \(t=0\) in the previous equation to obtain: \[ S(t=0)=S_0=e^{C_1}\cdot e^{-\lambda\cdot 0}=e^{C_1}\, \] ⊕The zero power of any number equals to one (\(e^0=1\)) so that our final formula: \[ S(t) = S_0e^{-\lambda t}\,. \]

To find how \(I(t)\) changes in time, we may use two ways. The first one, the simplest, is to notice that \(S(t)+I(t)\) is constant in time.Indeed, by using the dynamic equations, we get: \[ \frac{\mathrm d(S+I)}{\mathrm dt}=\frac{\mathrm dS}{\mathrm dt}+\frac{\mathrm dI}{\mathrm dt}=-\lambda S+\lambda S=0\,. \] It also means that \(S(t)+I(t)=S_0\). Then we obtain: \[ I(t)=S_0-S(t)=S_0-S_0e^{-\lambda t}=S_0(1-e^{-\lambda t})\,. \]

Another way is simply to substitute \(S(t)=S_0e^{-\lambda t}\) into the second dynamical equation for \(I(t)\). We get in this case: \[ \frac{\mathrm dI}{\mathrm dt} = \lambda S_0e^{-\lambda t}\quad\Rightarrow\quad \mathrm dI = \lambda S_0e^{-\lambda t}\mathrm dt=-S_0\mathrm de^{-\lambda t}\,. \] ⊕We check: \(\mathrm de^{-\lambda t}=\frac{\mathrm d}{\mathrm dt}(e^{-\lambda t})\mathrm dt=-\lambda e^{-\lambda t}\mathrm dt\). Integrating both sides: \[ \int\limits_0^t\mathrm dI =-S_0\int\limits_0^t\mathrm de^{-\lambda t}\quad\rightarrow I(t)-I(0) = -S_0(e^{-\lambda t}-e^{-\lambda\cdot 0})\,. \] Since \(I(0)=0\), we finally obtain: \[ I(t)=S_0(1-e^{-\lambda t}) \]

library(dplyr)
library(magrittr)

# Initial number of susceptibles: 100 people
S0 = 100

# Suppose that a person becomes sick after 2 weeks on average
lambda = 1/2

# We set a time to a sequence of one hundred points between 0 and 7 weeks
# This means that we observe a population for one week 100 times. 
# This will be our variable "t" in the data frame "df"
df = data.frame(t = seq(0,7,length.out=100))


# Then we add what we observe the following dynamics
df %<>% mutate(
  S = S0*exp(-lambda*t),
  I = S0*(1-exp(-lambda*t)))

# Library to show the tables nicely
library(knitr)

df %>%
  # Show first eight elements
  head(8) %>%
  # Show it nicely as a table
  kable(row.names=NA,digits=2,
    col.names=c("time, weeks","susceptibles",
                "infected"))

library(ggplot2)

df %>% ggplot(aes(x=t)) +
  geom_path(aes(y=I,color="Infected"),size=1.5) +
  geom_path(aes(y=S,color="Susceptibles"),size=1.5) +
  coord_cartesian(ylim=c(0,100)) +
  theme_bw(base_size=14,base_family='Times') +
  labs(x="Time (weeks)",y="Number of people") + 
  guides(color=guide_legend(title=NULL)) +
  theme(legend.position = c(1,.5),
        legend.justification = c(1, .5),
        legend.background = element_blank())

Example 2: Catalytic model with mortality

Now we add the mortality term into the \(I\)-compartment. We write: \[\frac{\mathrm dS(t)}{\mathrm dt} = -\lambda S(t),\quad \frac{\mathrm dI(t)}{\mathrm dt} = \lambda S(t)-\mu I(t),\] with the same initial conditions: \(S(0)=S_0\) and \(I(0)=0\).

After solving the first equation: \(S(t)=S_0e^{-\lambda t}\), we get the second equation in the form: \[ \frac{\mathrm dI}{\mathrm dt} = -\mu I+\lambda S_0e^{-\lambda t}\,. \]

Step 1. First we deal with the uniform solution. We solve: \[ \frac{\mathrm dI}{\mathrm dt} = -\mu I\quad\Rightarrow\quad \frac{\mathrm dI}I = -\mu\mathrm dt\quad\Rightarrow \mathrm d\ln I=-\mu\cdot\mathrm dt\,. \] If we integrate both parts of the equation, we get the solution dependent on the constant of integration \(C\): \[ I_{\text{uniform}}(t) = Ce^{-\mu t}\,. \]

Step 2. Here we want to obtain any particular solution of the original equation. To do so, we assume that the constant \(C\) varies in time: \(C=C(t)\), and we substitute this into the initial equation for \(I\): \[ \begin{split} \frac{\mathrm dI}{\mathrm dt} = \frac{d(C(t)e^{-\mu t})}{\mathrm dt} = \frac{\mathrm dC}{\mathrm dt}e^{-\mu t}+C\frac{\mathrm de^{-\mu t}}{\mathrm dt} = {} \qquad \\ \qquad \frac{\mathrm dC}{\mathrm dt}e^{-\mu t}-\mu Ce^{-\mu t}={-}\mu Ce^{-\mu t}+\lambda S_0e^{-\lambda t} \end{split} \] as we see the term \({-}\mu Ce^{-\mu t}\) is the same in both sides and it disappears. Then, we get: \[ \begin{split} \frac{\mathrm dC}{\mathrm dt}e^{-\mu t}=\lambda S_0e^{-\lambda t}\quad\Rightarrow\quad \frac{\mathrm dC}{\mathrm dt}=\lambda S_0e^{(\mu-\lambda)t} \qquad \\ \Rightarrow\quad \mathrm dC = \lambda S_0e^{(\mu-\lambda)t}\mathrm dt=\frac{\lambda S_0}{\mu-\lambda}\cdot\mathrm de^{(\mu-\lambda)t} \end{split} \] ⊕Again we check: \[ \frac{\mathrm de^{(\mu-\lambda)t}}{\mu-\lambda}=e^{(\mu-\lambda)t}\mathrm dt \]

By integrating both sides, we obtain: \[ C = \frac{\lambda S_0}{\mu-\lambda}e^{(\mu-\lambda)t} \] since we are interested in any particular solution, we can admit the constant of integration to any value22 Here we took: \(C(0)=1\). Thus, we write a particular solution in the form: \[ I_{\text{particular}}(t) = Ce^{-\mu t} = \frac{\lambda S_0}{\mu-\lambda}e^{(\mu-\lambda)t}\cdot e^{-\mu t}=\frac{\lambda S_0}{\mu-\lambda}e^{-\lambda t}\,. \]

Step 3.

Here we combine the uniform and particular into the final result: \[ I(t) = I_{\text{uniform}}(t)+I_{\text{particular}} = Ce^{-\mu t}+\frac{\lambda S_0}{\mu-\lambda}e^{-\lambda t}\,. \]

We remember that the initial condition \(I(0)=0\) should still work. So we use it to identify the constant \(C\): \[ I(t=0)=0=Ce^{-\mu\cdot0}+\frac{\lambda S_0}{\mu-\lambda}e^{-\lambda\cdot0}=C+\frac{\lambda S_0}{\mu-\lambda}\,, \] which means that: \[ C = -\frac{\lambda S_0}{\mu-\lambda}\,. \]

The final solution writes as follows: \[ I(t) = -\frac{\lambda S_0}{\mu-\lambda}e^{-\mu t}+\frac{\lambda S_0}{\mu-\lambda}e^{-\lambda t}\,. \] We can also re-write it in a bit nicer form: \[ I(t) = \frac{\lambda S_0}{\mu-\lambda}\left(e^{-\lambda t}-e^{-\mu t}\right). \] ⊕We may check that our solution gives the previous one if we take \(\mu=0\). Indeed: \[ I(t)\Big|_{\mu=0}=\frac{\lambda S_0}{0-\lambda}(e^{-\lambda t}-e^{-0\cdot t})=S_0(1-e^{-\lambda t}) \]

Notice that the number of died (or recovered) people \(R(t)\) will be simply expressed as follows: \[ R(t) = S(0)-(S(t)+I(t)) \]

# We add the mortality (or recovery) factor:
# each person recovers after 4 weeks on average
mu = 1/4

# We change our data frame a bit
df %<>% mutate(
  I = lambda*S0/(mu-lambda)*(exp(-lambda*t)-exp(-mu*t)))

# We also track the number of died (or recovered) persons
df %<>% mutate(
  R = S0-I-S
)

df %>% 
  head(8) %>%
  # Show it nicely as a table
  kable(row.names=NA,digits=2,
    col.names=c("time, weeks","susceptibles",
                "infected","recovered"))

df %>% ggplot(aes(x=t)) +
  geom_path(aes(y=R,color="Recovered"),size=1.5) +
  geom_path(aes(y=I,color="Infected"),size=1.5) +
  geom_path(aes(y=S,color="Susceptibles"),size=1.5) +
  coord_cartesian(ylim=c(0,100)) +
  theme_bw(base_size=14,base_family='Times') +
  labs(x="Time (weeks)",y="Number of people") + 
  guides(color=guide_legend(title=NULL)) +
  theme(legend.position = c(.5,1),
        legend.justification = c(.5, 1),
        legend.background = element_blank())

df %>% 
  # Get the last row of the data frame
  tail(1) %>% 
  # Show it nicely as a table
  kable(row.names=NA,digits=0,
    col.names=c("time, weeks","susceptibles",
                "infected","recovered"))