Ryan Johnston

Pandemics are a widespread, rapid spread of disease, with exponentially rising cases over a large area. Endemic viruses, meanwhile, are constantly present and have a fairly predictable spread. That predictability allows health care systems and doctors to prepare and adapt, reducing loss of life.

The coronavirus (COVID-19) pandemic was one of the most globally disruptive events in decades. It first appeared in Wuhan, China during December of 2019 and quickly spread across the entire world in a matter of months. As of February 9th, 2023, there are approximately 672.448 million confirmed cases and 6.85 million deaths globally according to John Hopkins University. After two years, the World Health Organization says Covid-19 still remains a global health emergency, but the pandemic is now at a transition point, possibly towards an endemic phase. The affects of the the pandemic are still being felt, with supply chains and economies struggling to recover and healthcare systems at their breaking point.

There has been an extensive amount of research that has been done over the past year on COVID-19, including the origin and cause, the spread of the virus, treatment, and a possible cure or at the very least a vaccination. Additionally, there has also been a large number of studies investigating non-medical aspects related to COVID-19 such as it’s impact on the economy, the role of politics, religion, and how it has changed human culture. However, there are still plenty of unknowns that remain as the COVID-19 pandemic continues into a possible endemic phase in 2023 and 2024. There have now been over 13.28 billion doses of vaccines administered to protect individuals from the adverse of COVID-19 and long-COVID.

In this project we use random-walk and Monte-Carlo methods to simulate the spread of virus. We explore four different scenarios that account for different parameters such as lockdowns and varying mobility. We then improve our code by accounting for hospitalizations and the vulnerable section of the population. Later on, we explore the spread of a fictional virus with some interesting results.

Overview

We used Random Walk methods in order to study the stochastic behavior of the spatio-temporal propagation of a virus. The model consists of individuals that can interact with their neighbors within a range and at a certain mobility, with a certain probability of interaction, which in turn can produce infection with the result of recovery or death. The approach relies on physics concepts similar to those you used for the radiative transfer problem. We compare the mortality rate from your simulations to that of a normal flu (i.e \(\sim 0.1\%\)) and covid-19 ( \( \sim 1\%\)).

We start with a two-dimensional square lattice or grid of N by N squares with one individual per square; indices are \(i = 1,N\) and \(j = 1,N\). We will assign 4 possible states (each taking up a value of 0 or 1) to each individual as the virus spreads:

1. Susceptible (to infection)
2. Infected
3. Recovered
4. Dead

Once an individual located at coordinate \(k = (i, j) \) is infected it can affect other individuals around (i.e. spread the virus) and the process continues until, in principle, an endemic or a pandemic is reached1. We start the simulations by initially assuming that the central cell is occupied by an infected individual; i.e. ( \( i(N=2,N=2) = 1\) ) while all of the other individuals are susceptible (i.e. \( s(k) = 1\) while \(i(k) = r(k) = d(k) = 0 \)).

Inverse Sampling

The fundamental principle is defined to be \[ \psi(x_0) = \frac{\int_{a}^{x_0} P(x) dx}{\int_{a}^{b} P(x)dx} \] which in the case that \(1 = \int_{a}^{b} P(x)dx\) can be re-written as \[ \boxed{\xi = \int_{a}^{x_0} P(x) dx = \psi(x_0)}\] where \(\xi\) is a random number sampled from a uniform distribution from 0 to 1. We can apply the fundamental principle to find the inverse cumulative distribution functions for two key distributions.

The algorithm

The basic idea behind the algorithm is to first estimate what is called the range (or zone) of influence for each individual. Then within this zone we estimate how many individuals the individual at \(k\) interacted with. These are the susceptible individuals. We consider that each one of them can be infected with a probability \(0 \leq p_{i} \leq 1\).

We start at t = 0 and take 1 day time-steps. We then apply the following for each cell containing an infected individual

Zone of influence probability

The zone of influence probability, \(p_{zi}\): Cell \(k\) can infect only people in other cells within a finite region of size \(R_{zi}\). This is the zone of influence in a given day. Now the simulation starts by computing the value of Rzi randomly with some probability distribution. We assume an exponential law used in the literature

Now the simulation starts by computing the value of \(R_{zi}\) randomly with some probability distribution. We assume an exponential law used in the literature. So, the zone of influence is defined to be \[ p_{zi} = \frac{e^{-R_{zi}/R_{0}}}{R_{0}} \tag{1}. \] which indicates that is less probable to move far away from the average interaction range \(R_{0}\) of the individuals. To simplify the algorithm in Cartesian coordinates we assume that the region is a square of side \(2R_{zi}\) (the closest integer to \(2R_{zi}\)), centered at the cell \(k\). This is one of the parameters of the model. First, we check if the probability density function (PDF) is normalized. \[ \begin{align*} 1 &= \int_{0}^{\infty} \frac{e^{-R_{zi}/R_{0}}}{R_{0}} dR_{zi} \\ 1 &= \left[-e^{-R_{zi}/R_{0}} \right]_{0}^{\infty} \\ 1 &= 1 %-\frac{A^{2}}{2}e^{-2R_{zi}/R_{0}} \end{align*} \] So, Equation (1) is indeed normalized. Next we compute the cumulative distribution function (CDF) as follows \[ \begin{align*} \xi\left(R_{zi}\right) &= \int_{0}^{R_{zi}} \frac{e^{-R_{zi}/R_{0}}}{R_{0}} \; dR_{zi} \\ &= \left[-e^{-R_{zi}/R_{0}} \right]_{0}^{R_{zi}} \\ &= 1 - e^{-R_{zi}/R_{0}} \end{align*} \] Solving for \(R_{zi}\) gives \[ \begin{align} \boxed{R_{zi} = -R_{0}\log\left[1-\xi\right]} \end{align} \]

Lockdown

A lockdown is simulated by lowering the value of \(R_{0}\). In other words, reducing the zone of influence. We assume that lockdown starts at \(t_{lock} = 3\) weeks, or 21 days, since the start of the simulation at \(t = 0\).

Mobility Probability

\(p_{m}\): Once \(R_{zi}\) is chosen for an infected individual, it can interact randomly with any of the individuals within the zone of influence. We characterize the interaction between individuals by a cross section \(\sigma\); it measures the probability of interaction between two individuals during a given day. The number of possible interactions is then given by the probability formula \[ N_{int} = m \times n\sigma \] where \(n\) is the surface density (naturally; n = 1 our model). The parameter m defines the mobility of the individual within its zone of influence and it can be thought of as the temperature of the system. This allows us to use statistical physics and assume the mobility distribution is of the form \[ p_{m}(m) = P(n) = \frac{e^{-m/T}}{e^{1/T}-1} \] where the parameter \(T\) is the temperature of the system; the hotter the system the higher the mobility. In reality, the optical-depth is \(\tau = n\sigma (2R_{zi})\) for an individual moving a distance \(2R_{zi}\) in a medium of volume density \(n\) (\(cm^{3}\)) with an interaction cross-section \(\sigma\) (\(cm^{2}\)). However, because we are using a 2-dimensional grid and the density is per unit area (\(cm^{2}\)) and not per unit volume, it makes more sense to talk about a dimensionless mobility parameter rather than speed in order to keep \(\tau \) dimensionless. Now, we must normalize the probability distribution function \begin{align*} 1 &= A \int_{0}^{\infty} p_{m}(m) dm \\ 1 &= A \int_{0}^{\infty} \frac{e^{-m/T}}{e^{1/T}-1} dm \\ 1 &= A \int \left[-\frac{T e^{-m/T}}{e^{1/T}-1} \right]_{0}^{\infty} \\ 1 &= A \left(\frac{T}{e^{1/T}-1}\right) \end{align*} Solving for the normalization constant \(A\) gives \[ A = \frac{e^{1/T}-1}{T} \] So, the normalized PDF is \[ \boxed{p_{m}(m) = \frac{e^{-m/T}}{T}} \] which is of a similar form to that of Equation (1). The CDF is then computed as \[ \begin{align*} \xi\left(m\right) &= \int_{0}^{m} \frac{e^{-m/T}}{T} dm \\ &= \left[-e^{-m/T} \right]_{0}^{m} \\ &= 1 - e^{-m/T} \end{align*} \] Solving for \(m\) then gives \[ \boxed{m = - T\log\left[1 - \xi\right]} %\left({m}\right) \]

The infection probability

pi: Once we know the range Rzi and the number of interactions \(N_{int}\) of an infected cell, we choose randomly \(N_{int}\) cells, with coordinates \(k_{0}\), within the interaction region. In other words we randomly chose \(N_{int}\) susceptible individuals and we:

1. Assign to them a value \(s(k') = 1\).
2. Decide randomly if it is infected or not with with probability \(0 \leq pi \leq 1\).
3. The daily increments, \(\Delta S(t)\), \(\Delta I(t)\), \(\Delta R(t)\), \(\Delta D(t)\).
4. Do two things in case that it becomes infected:
• We change the value of the corresponding matrix elements \(s(k')= 0\) and \(i(k') = 1\).
• We store the time it becomes infected in the infection time matrix \(TM_{i}(k) = t\). This matrix is important for the next step in deciding when the individual becomes recovered or dead.

The death probability

\(p_{d}\): After interacting with all the \(N_{int}\) individuals within the range \(R\), during the day, finally we decide if cell \(k\) becomes recovered (stays in the game) or dead (removed from the grid) at the end of the day.

We wait a time \(t_{r} = 21\) days (from the day of infection which you should have stored in the infection time matrix \(TM_{i}(k))\) and we decide randomly with probability \(p_{d} \leq 1\) if \(k\) dies from the infection. Finally we store the matrices \(r(k) = 1 (or 0)\), and \(d(k) = 0\) (or 1) accordingly. In case \(k\) is dead, we set that cell to \(i(k) = 0\).

Scenario Parameters

There are 7 parameters in total: \( N \), \(R_{0}\), \( \sigma \), \(T\), \(t_{r}\), \(p_{i}\) and \(p_{d}\). We apply our algorithm to the cases shown in Table 1 in order to portray the effect of lockdown (\(R_{0}\) effect) and mobility (\(T\) effect) on the spread of the virus. You will be computing the state of the system over 100 days with a time-step of \(\Delta t = 1\) day. At the end of each day you should store:

1. In the matrices \(s\), \(i\), \(r\), \(d\) the current state of the system.
2. The total (or cumulative) number of susceptible \(S(t)\), infected, \(I(t)\), recovered, \(R(t)\) and dead, \(D(t)\)
3. The daily increments, \(\Delta S(t)\), \(\Delta I(t)\), \(\Delta R(t)\), \(\Delta D(t)\).

Lockdowns kick in at \(t_{lock} = 21\) days (3 weeks) where \(R_0\) is then reduced from 8 to 3.

Our first task is to explore the cumulative death counts and daily death rates over time. This is shown in Figure 1 in a sim- ple case with no lockdowns for 200 days. For the daily death counts we also compute the rolling 7-day average. We can see in Figure 1 that the daily death rate increases dramatically and then plateaus after about 80 days. The cumulative death rate increases linearly after this time as well. This is indicative of a transition from pandemic to endemic, where the population is not longer affected by spread of the virus. In other words, the virus is regularly occurring within an area.

For each scenario, the spread of the virus was simulated for 115 days. Snapshots of the virus’ spread were generated every 14 days (Listing 18). We also created a few .MP4 animations that show the spread of the virus day by day (Listing 19). Figures 2, 3, 4, and 5. shows eight 14 week snap shots of the spread of the virus for Scenarios 1 through 4, each for a duration of 115 days. In scenarios 2 and 4 the lockdowns that are implemented after 21 days is very effective at slowing the spread of the virus, and would re- duce the burden on the healthcare system. Despite having no lockdowns in Scenarios 1 and 3, the lower mobility in Scenario 3 makes the virus spread much slower.