Table of Contents

• Abstract
• 1. Introduction
• 2. Models and Results
  • 2.1 Basic SIR
  • 2.2 Spatial SIR
    • 2.2.1 Agent-based Model
    • 2.2.2 PDE Model
  • 2.3 SIR with varition on Parameters
    • 2.3.1 Agent-based Model
    • 2.3.2 ODE Model
  • 2.4 SEIR
    • 2.4.1 Agent-based Model
    • 2.4.2 ODE Model
  • 2.5 Fit to Data
• 3. Discussion
  • 3.1 Conclusion
  • 3.2 Computational Performance
  • 3.3 Limitation
  • 3.4 Further Investigation
• 4. Bibliography

Abstract

In this report, we study the spread of a disease using the SIR model and its variations. We start with the basic SIR model to help us choose interesting parameters. We then explore the spatial SIR model with the set of chosen parameters, where we solve the model using the agent-based model and the PDE model, and study the effect of the diffusion parameter and the starting position on the evolvement of the infectious population. We then move on to the three variations we proposed in the midterm checkpoint. Finally, we present our findings by studying these models and our discussion on the limitations and improvements of these models.

1. Introduction

To model the spread of a new disease, we classify individuals into three categories according to the SIR model:

1. Susceptible - individuals who have not yet caught the disease.
2. Infectious - individuals who are sick and may spread the disease to susceptible individuals.
3. Recovered - (aka. removed) individuals who were previously infectious, and either have recovered and are now immune, or have died. Either way they can not get the disease again or infect susceptible individuals.

We use $S(t)$, $I(t)$ and $R(t)$ to denote the number of individuals in each category above at time $t$, and $N$ to denote the total number of individuals, so $S(t) + I(t) + R(t) = N$ at all times. We define the corresponding fractions for each group by $s(t) = S(t) / N$, $i(t) = I(t) / N$, and $r(t) = R(t) / N$, so $s(t) + i(t) + r(t) = 1$ at all times.

In the plots blow, we will use blue for susceptible, orange for infectious and green for recovered.

2. Models and Results

2.1 Basic SIR

We use the Kermack-McKendrick Model which has two parameters to model the change in $S(t)$, $I(t)$ and $R(t)$:

• $b$: the number of contacts that an infectious individual has each day
• $k$: the fraction of the infectious population which recovers each day

At time $t=0$, there is a small fraction of infectious population $i_0$. Each day an infectious individual randomly contacts $b$ other individuals, who can be from one of the three categories but only susceptible individuals get infected. Meanwhile, $k$ fraction of the infectious population recovers.

We use two methods to find how $s(t)$, $i(t)$ and $r(t)$ change across time with different values of $b$ and $k$:

1. agent-based model, where we represent a person using a class with an internal state which is one of $S$, $I$, or $R$.
2. ODE model, where we model the fraction of each population through the following system of differential equations: \begin{aligned} \frac{d s}{d t} &= -b\cdot s(t)\cdot i(t) \\ \frac{d i}{d t} &= b\cdot s(t)\cdot i(t) - k\cdot i(t) \\ \frac{d r}{d t} &= k\cdot i(t) \end{aligned}

We get the following phase diagram of the infectious rate $i(t)$ by $b,k$ using the agent-based model in the midterm checkpoint:

In particular, we can take a look at three timestamps during the spread of the disease

The results of the ODE model are similar but the colors change more smoothly. This is because unlike the agent-based model, there is randomness in the ODE model.

We now choose $b$ and $k$ near the phase transition boundary for in-depth analysis in the later sections.

From the plots we can see that for the late stage such as $t = 50$, when $b$ does not matter that much anymore, $k = 0.05$ seems to be the boundary of completely white and a little colorful, so we will choose this $k$. For the early stage such as $t = 5$, when $b$ does have an effect, the boundary seems to be around 5, so we choose $b = 7$ for the agent-based model and $b = 5$ for the PDE model.

2.2 Spatial SIR

2.2.1 Agent-based Model

We model the spatial SIR by adding a 2D position coordinate $(x,y)$, where $0\le x \le 1$ and $0 \le y \le 1$. At each day, each agents moves with length $p$ in a random direction, and infected people will infect others within radius $q$ . We use KDTree to find such $r$ nearest neighbors for each infected people. For simulation, we set $N = 10000, i_0 = 0.1, b = 7, k = 0.05$, which results $q = 0.015$.

Effect of step size $p$

To investigate how infected will change over $p$, we run 20 simulations with $p$ equally spaced between 0 to 1. The results of $s(t)$, $i(t)$ and $r(t)$ change over $p$ is shown below.

The depth of color indicates the magnitude of $p$ The darker color indicates larger value $p$ From the plot, we can see that with $p$ starting at 0 (the lightest line,) the $s(t)$ lines skew toward right as $p$ increases, and the peak of the infectious rate also increases. The max of the peak reaches at approximately $p = 0.5$, with $i$ approximately equal to 0.8. Then the lines start skew toward left, and the peaks start to decrease again. This result is due to we limit the agent mobility by setting them stay at the same position if their move is out of bound. If $p$ is too large, agents are equivalent to not move. Thus both of the spread and the maximum of infectious rate will decrease.

It is also interesting to note that the susceptible fraction $s(t)$ and the recovered fraction $r(t)$ change the fastest with certain $p$ around 0.5, and change the slowest for too small and too large $p$'s. As a result, we suspect that there exists an optimal $p$ such that the agents in all three categories moves the fastest on average and hence all three fractions $s(t), i(t), r(t)$ change the fastest. To find such $p$, we can find the largest averaged actual step size $\bar p$ for all agents in a population (recall that some agents stay in their position if they would step outside the domain).

More specifically, for each value of $p$ from 0 to 1.5, we

• simulate 100 populations, each with 500,000 agents
• for each population
  • use the method move(p) to obtain the updated agent positions
  • find the difference in their positions, which is the actual step size
  • take average for all agents in the population
• take average for all 100 populations

The averaged actual step size changes with $p$ in the below fashion

Thus, the optimal $p$ should be around 0.49, with which the agents move the fastest and all three fractions change the fastest.

Effect of starting position

To investigate how starting position of infected will affect the spread, we choose three $p = 0, 0.5, 1$. For each $p$ we run total number of 30 simulations with 5 sets starting at center, 5 sets starting at corner and 20 sets starting random. The resulted plots are shown below, with $p =0, 0.5, 1$ from left to right.

By comparing the plots, we can conclude that if $p$ is within the middle range, ie. 0.3 - 0.5, the starting positions really does not matter, since the spread of infectious is so quick. The middle point $p$'s indicate agents can freely move around in any direction. Thus, even starting at the corner, the infected people are likely to get to the center at anytime point. With $p$ really small (we set extreme to 0 in this case,) since corner infected can only infect 1/4 of the people around, the spread is much lower than center/random starts. What's more, we can see with some extreme cases, infected might even not be able to spread the virus, with all of them getting recovered without further infection. We also notice that no matter what starting position is, lower mobility (either extreme large $p$ or small $p$) results lower infections rate and spread. Thus, stay at home is indeed important to flatten the curve.

To verify our conclusion, we further visualize the spread of infectious by scatter plots. The below shows an extreme case that with $p = 0$, starting at corner can hardly spread the disease if an infectious agent has no susceptible neighbors in radius $q$. The recover rate outspeed the spread rate. While starting at center and random, it is more likely to spread the disease, since there are more contacts within the radius.

With $p = 0.5$, agents move the fastest and all most no people are susceptible within 10 days regardless of the starting points.

When $p = 1$, the spread is slower than the case when $p = 0.5$. Yet, the peripheral agents still have mobility.

2.2.2 PDE Model

The PDEs we want to solve are: \begin{align} \frac{\partial s(\mathbf{x},t)}{\partial t} &= -b\cdot s(\mathbf{x},t) \cdot i(\mathbf{x},t) + p\cdot \Delta s(\mathbf{x},t) \\ \frac{\partial r(\mathbf{x},t)}{\partial t} &= k\cdot i(\mathbf{x},t) + p\cdot \Delta r(\mathbf{x},t) \\ \frac{\partial i(\mathbf{x},t)}{\partial t} &= b\cdot s(\mathbf{x},t)\cdot i(\mathbf{x},t) - k\cdot i(\mathbf{x},t) + p\cdot \Delta i(\mathbf{x},t) \end{align}

We want to use finite differencing to transform these PDEs into systems of ODEs and then use solve_ivp in the scipy package to get the solutions. We first discretize $[0,1]$ on the $x$-axis to $0 = x_0 < x_1 < \cdots < x_M = 1$ where $x_i - x_{i-1} = \frac{1}{M} := h$ for all $i$ and similarly for $y_0,\ldots,y_M$.

Consider the first PDE on $s(\mathbf{x},t)$. The Laplacian at $(x_j,y_k)$ becomes \begin{align*} \Delta s(x_j,y_k,t) &= \frac{\partial^2 s(x_j,y_k,t)}{\partial x^2} + \frac{\partial^2 s(x_j,y_k,t)}{\partial y^2} \\ &= \frac{s(x_{j+1},y_k,t) - 2s(x_j,y_k,t) + s(x_{j-1},y_k,t)}{h^2} + \frac{s(x_j,y_{k+1},t) - 2s(x_j,y_k,t) + s(x_j,y_{k-1},t)}{h^2} \end{align*}

Define $s_{j,k}(t) = s(x_j,y_k,t)$. Then the PDE becomes a system of ODEs: \begin{align} \frac{ds_{j,k}(t)}{dt} = -b\cdot s_{j,k}(t)\cdot i_{j,k}(t) + \frac{p}{h^2} [s_{j+1,k}(t) + s_{j-1,k}(t) + s_{j,k+1}(t) + s_{j,k-1}(t) - 4s_{j,k}(t)], \quad \forall j,k = 1,\ldots,M \end{align}

Similarly, we define $i_{j,k}(t)$ and $r_{j,k}(t)$ and the other two PDEs become \begin{align} \frac{dr_{j,k}(t)}{dt} &= k\cdot i_{j,k}(t) + \frac{p}{h^2} [r_{j+1,k}(t) + r_{j-1,k}(t) + r_{j,k+1}(t) + r_{j,k-1}(t) - 4r_{j,k}(t)], \quad \forall j,k = 1,\ldots,M \\ \frac{di_{j,k}(t)}{dt} &= b\cdot s_{j,k}(t)\cdot i_{j,k}(t) - k\cdot i_{j,k}(t) + \frac{p}{h^2} [i_{j+1,k}(t) + i_{j-1,k}(t) + i_{j,k+1}(t) + i_{j,k-1}(t) - 4i_{j,k}(t)], \quad \forall j,k = 1,\ldots,M \end{align}

To write these in matrix form, we need to define the second-order finite differencing matrix \begin{align*} D = \begin{bmatrix} -2 & 2 & 0 \\ 1 & -2 & 1 & 0\\ 0 & 1 & -2 & 1 & 0 \\ & & \ddots & \ddots & \ddots \\ & & & 0 & 1 & -2 & 1 \\ & & & & 0 & 2 & -2 \end{bmatrix} \end{align*} where the empty entries are all 0. Then we can write the three system of ODEs as \begin{align*} \frac{ds(t)}{dt} &= -b\cdot s(t)\times i(t) + \frac{p}{h^2}(D\cdot s(t) + s(t)\cdot D^\top) \\ \frac{di(t)}{dt} &= k\cdot i(t) + \frac{p}{h^2}(D\cdot i(t) + i(t)\cdot D^\top) \\ \frac{dr(t)}{dt} &= b\cdot s(t)\times i(t) - k\cdot i(t) + \frac{p}{h^2}(D\cdot r(t) + r(t)\cdot D^\top) \end{align*} where $s(t),i(t),r(t)$ are the $M\times M$ matrix-valued functions and $s(t)\times i(t)$ denotes the entry-wise multiplication and $D\cdot s(t)$, etc., denote the matrix multiplication.

How do simulations change with $p$?

We choose, according to the experiments in the midterm checkpoint, $M = 200$, $i_0 = 0.001$, $T = 100$, $b = 5$, and $k = 0.05$.

We run the simulation with $p$ in np.logspace(-4, -1, 10) / M and draw the line plots for the mean of $s(t), i(t), r(t)$ where the mean is taken over all the grids. The plots look like (where darker lines correspond to large $p$'s):

We can see that $p$ increases, the number of infectious people reaches its peak earlier and faster. This makes sense because in the PDEs, $p$ is the weight of the Laplacian term, i.e., the spatial diffusion. Hence we should expect that larger $p$ results in faster spread of the disease.

Now since the rate of spread is uniformly increasing with $p$, we will choose a moderate $p = 0.001 / M$ so that the simulation can be finished in reasonably short time and the timestamp plots can render a clear evolution of the different populations.

How do the three populations, with different starting points, evolve with time?

We fix $p = 0.001 / M$.

If we start from the center:

If we start from the corner:

If we start from a random point:

We can see that at first, the infectious population spreads out from the starting point and the susceptible population fades away against the infectious population; then, the recovered population starts to increase, again from the starting point, and chases the infectious population. In the end, all population would become recovered.

When the starting point is the center, or some random points near the center, the rate of spread is very high such that by time $t = 50$, most of the population are first infected and then recovered. However, when the starting point is the corner, the rate of spread decreases so that by time $t = 50$, only about a quarter of the population are infected and then recovered.

How does the infectious group change with different starting points?

We can see that starting from the center gives much faster spread than the corner, and starting at random points sometimes gives, when the random point is near the center, similar results as starting from the center, and sometimes gives, when the random point is near the corner or the side, similar results as starting from the corner.

2.3 SIR with varition on Parameters

The basic SIR model has some unrealistic settings.

• The interaction rate $b$ can also change over time. For example, $b$ decreases with the enforcement of social distancing and quanrantine and increases (mostly locally) with mass gatherings. The accumulated effect of different actions can be hard to measure. Therefore, we can model $b$ using a piece-wise constant function such as \begin{align*} b(t) = \begin{cases} b_1, & t < t_0 \\ b_2, & t_0 \le t \le t_1 \\ b_3, & t \ge t_1 \end{cases} \end{align*} where $t_0$ represents the time that hospital treatments of the disease become mature, and $t_1$ represents the time that the vaccine is invented, and $b_1<b_2<b_3$.

• Similarly, the recover fraction $k$ may depend on time. As recovery methods, such as vaccine and hospitalization, are put into practice, the rate of recovery is likely to increase. so we can again use a piece-wise constant function to model $k$.

• In addition, we assumed that an infectious individual will spread the disease to all the susceptible people they meets, but generally that's not the case, since the infection can depend on how long they stay together, whether they wear masks, the air ventilation around them, etc. Therefore, we can add an additional parameter $p$ to model the probability that an infectious individual actually spreads the disease when he has contact with an susceptible individual.

We will investigate how $s(t)$, $i(t)$ and $r(t)$ change with these time-dependent parameters, using an agent-based model and an ODE model. In particular, we consider the following three effects:

• Social Distancing

  Social distancing can affect $b$, the number of interactions made by an infectious individual per day. Since social distancing rules take time to be implemented, we use

\begin{align*} b(t) = \begin{cases} 10, & t \le 5 \\ 3, & 5 < t \le 30 \\ 1, & 30 \le t \le 100 \end{cases} \end{align*}

• Drug Development and Distribution

  This will affect the removed fraction per day $k$. We use $k=0.01$ when there is no drug, and $k=0.1$ for the early distribution period of the drug, and $k=0.8$ for the distributed period.

\begin{align*} k(t) = \begin{cases} 0.01, & t \le 15 \\ 0.1, & 15 < t \le 70 \\ 0.8, & 70 \le t \le 100 \end{cases} \end{align*}

• Virus Mutation

  If mutation happens, the virus becomes more infectious. This can be measured by the transmission probability $p$. Holding other constant, if two individuals has an interaction, then we expected $p$ to be higher if the virus becomes more infectious. In the experiments below we use $p=0.1$ for the normal virus first 20 days and $p=0.8$ for the mutated virus in the last 80 days.

\begin{align*} p(t) = \begin{cases} 0.1, & t \le 20 \\ 0.8, & 20 \le t \le 100 \end{cases} \end{align*}

2.3.1 Agent-based Model

We first make a plot as a base case to be compared to, where $b=10, k=0.01, p=0.1$ We see from the plot below that $i(t)$ reaches its maximum at around 80% on day 30.

Let's see the effect of social distancing. The trajectories of $s(t), i(t)$ and $r(t)$ for altered parameters are represented in dash lines. We see that early social distancing (left) can mitigate the spreads of the pandemic. The dash orange line $i(t)$ increases at a slower speed. The peak arrives later with a smaller value, which could alleviate the pressure of medical capacity. In contrast, late social distancing starting at day 20 (right) doesn't really help, since the susceptible fraction $s(t)$ is already quite small and most people are infected or removed.

Then, we analyze the effect of drug development and distribution. Similarly, early drug distribution from day 15 (left) effectively change the trajectory of $i(t)$: it starts to decreases immediately rather than continue to reach its maximum in the base case. However, the susceptible fraction $s(t)$ does not change much. Drug distribution in a late stage around day 30 (right) changes $i(t)$ immediately too, but the peak is already arrived and the medical system may already be crashed.

For the effect of virus mutation, we see that if the mutation happens in a pre-peak stage, then the peak will arrive earlier with a higher value. In contrast, if the mutation happens in a post-peak stage around day 80, then nothing changes since the susceptible fraction is already 0 – no more individual can get infected, and every individual is either infected or removed.

We can consider the three effects in the early stage simultaneously. The dash curves deviate significantly from the base case. Social distancing rules are enforced on day 5, drugs distribution starts from day 15. Theses two methods effectively decreases $i(t)$. The virus mutation happens on day 20, where we see a jump in the orange dash line. But the peak value is still smaller than that in the base case. Unfortunately, all susceptible individuals are infected eventually.

2.3.2 ODE Model

We run similar simulation using the ODE model. One thing difficult is how to put the parameter $p$ in to the differential equations. Since $p$ is related to the speed of virus, transmission, we can just substitute $b$ in the basic SIR model by $b\times p$.

But we found that with the base case parameter $b=10, k=0.01, p=0.1$ the ODE model generates a different base case plot from the agent model's. We may consider adding a power to $p$ to make the two base case plots similar. After tunning the power we let it be 1.5.

The new differential equations are

\begin{aligned} \frac{d s}{d t} &= -b\cdot p^{1.5} \cdot s(t)\cdot i(t) \\ \frac{d i}{d t} &= b\cdot p^{1.5}\cdot s(t)\cdot i(t) - k\cdot i(t) \\ \frac{d r}{d t} &= k\cdot i(t) \end{aligned}

The effects of the three time-dependent parameters are shown below.

We see that social distancing seems more significant in the ode model than in the agent model. All three fractions changes quite slowly. The mutation effect is more significant too with a large increment in $i(t)$ on the mutation day $t=20$. Maybe it is due to the effect of introducing the power term in $b\cdot p^1.5$. Drug distribution have similar effects in both models.

In all, from the experiments we conclude that, by affecting the interaction rate $b$, early social distancing can flatten the $i(t)$ curve, postpone the arrival of its peak with a smaller value. Drug distribution facilitates the process of recovery by changing $k$. These two methods affect $i(t)$ from different way: one makes the increment of $I$ from $S$ smaller, one makes the decrement of $I$ to $R$ larger. The mutation effect can significantly speed up the transmission if it happens in an early stage.

The takeaway is that, social distancing and drug distribution should start as early as possible to mitigate the spread of the pandemic.

2.4 SEIR

In reality, there is a latent period of Covid-19, when an individual has been infected but is not yet infectious. As a result, we may consider adding a new state named exposed. This is the SEIR model.

The system of differential equations becomes \begin{aligned} \frac{d s}{d t} &= -b\cdot s(t)\cdot i(t) \\ \frac{d e}{d t} &= b\cdot s(t)\cdot i(t)- f\cdot e(t) \\ \frac{d i}{d t} &=a\cdot e(t) -k \cdot i(t) \\ \frac{d r}{d t} &=k \cdot i(t) \end{aligned}

where $f$ is the fraction of the exposed individuals that becomes infectious each day.

We can repeat our simulation above to see how the curves of $s,i,r$ change with different values of $f$. And plot the 2-d and 3-d phase diagram of the three parameters $b, f, k$. Note that we need to specify the initial fraction of exposed population $e(0)$.

2.4.1 Agent-based Model

Following the coloring tradition, we use blues, oranges, greens for $S, I, R$. In addition, we use purples for the new state $E$.

We set $e(0) = i(0) = 0.001, b = 1, k = 0.01$ and run experiments for 10 different $f$ values in log-spaced interval 0.01 to 1. Deeper lines corresponds to larger values of $f$.

It can be derived that when $f$ is getting closer to 1, the SEIR model converges to the SIR model since the exposed period is quite short. This is shown in the top right purple plot of $e(t)$, where the curves becomes more flattened as $f$ increases.

When $f$ is very small, the individuals "stuck" in the state $E$ and transit to state $I$ very slowly. This is good since we flatten the orange $i(t)$ curves which gives the medical system more time to prepare and react.

We then investigate the phase diagrams of the parameters $b,f,k$. We simulate $i(t)$ for 100 days with 1,000 $(b,f,k)$ combinations. The 3-D phase diagram on the 50th day is shown below. Again, darker dots represents larger values of $i(50)$. After investigation we found that $i(50)$ is not sensitive to $k>0.3$, which is consistent with our findings in the SIR model. So we only show the results for small $k$.

From the above results we can generate three the 2-d phase diagrams by taking the mean of $i(t)$ in the remaining dimension.

We can view each cell as a population instance, for different combinations of parameters. For population instances with larger infected fraction $f$, the peak of $i(t)$ (darker orange) arrives earlier. If $f$ is small, then the peak arrive later since many individuals are stuck in the exposed state. The peak value also varies. The population instances with larger recovery rate $k$ tends to recover faster, and the color is lighter. The interaction rate $b$ has no significant effect when $t$ is large, which is consistent with the results in the midterm checkpoint.

2.4.2 ODE Model

Similarly, we repeat these experiments using the ODE model.

The results agree with those of the agent-based model.

2.5 Fit to Data

In this section we try to fit out ODE models to real life Covid-19 data.

The number of daily new cases can be retrieved from CDC's website.

Unfortunately, there is no data for daily removed cases which is necessary to estimate the removed rate $k$. For simplicity, we assume $k=0.05$, and use the number of daily new cases on day $t-1/0.05=t-20$ to estimate the number of daily removed cases on day $t$.

Therefore, from the data of daily new cases $N(t)$ we can compute $$S(t)=1-\sum_{s\le t}N(s)$$ and $$I(t)=\sum_{s\le t}[N(s)-N(s-20)]$$

Dividing these two quantities by the US population of 382 million we can get the plot of susceptible fraction and infectious fraction.

The two curves deviate from those in our basic ODE model in previous sections. This is because the population size is large and the parameters $b,k$ are affected by various factors, such as social distancing policies, as discussed in the section of SIR with variation on parameters.

As a result, we may consider a piece-wise constant function $b$. The cutoff points can be chosen from the days of local minimums and the local maximums, i.e. day 94, 142, 191, and 239. For each period we use solve_bvp and pass the boundary conditions. Note that we start from the 50-th day since the two lines are nearly constant and the ode result is be unstable in the earlier period.

We can see that the fitting to $s(t)$ is good, while there are rooms for improvement for fitting $i(t)$. The estimated interaction rate parameter $b$ is large from day 50 to 94, and decreases and increases again in the summer (day 150-200) and the election periods (day 250+). Also note that the values is quite smaller than those used in the simulation in the previous sections. But such small $b$ already has a huge impact on the society. We can also project that if $b$ is a constant from day 300, then the curve goes exponentially again.

For a better fitting, we can cut the curves into more pieces, especially near the local minimums/maximums, where $b$ should change quickly. We may also consider variation on the removed fraction $k$ if the relevant data is available.

3. Discussion

3.1 Conclusion

Regarding the two questions asked in spatial.md for each of the two models, we have the following conclusions:

• the role of $p$ differs for the agent-based model and the PDE model. In PDE, $p$ is the weight of the Laplacian term, so larger $p$ corresponds to faster spread. In agent-based model, however, the rate of spread would first increase and then decrease as $p$ increases. This is because too large a $p$, similar as too small a $p$, restricts the movement of the population. For example, if $p > \sqrt{2}$, then all infectious individuals would never move away from where they are.

• if we fix the interesting $p$ (again, the interesting $p$ is different for the two models), we would observe that the initial infection at the center is much more effective than the initial infection at the corner. This makes sense because infectious individuals at the center would have more space to move around and infect more people, while infectious individuals at the corner are kind of limited in their movements.

3.2 Computational Performance

In this section we discuss how we applied the computational thinking learned from this course into practice to optimize our algorithms and facilitate the simulation.

Agent-based model:

We notice that if we use the method of creating a Person() object for each individual, the simulation is quite slow. For instance, if the population size $N=10000$, we will have 10000 objects in the memory and need to loop over each of them. As a result, we continuously find ways to optimize it.

• in the final report, we create a Population() object instead to represent a population. The information to store is just three sets of indices for the three population $S, I$, and $R$. In simulation, we use methods pop.infect(b), pop.recover(k) to update the sets of indices. This turns out to be much 6x faster than the Person() method when $N=10000$.

• in order to simulate the interaction, we use np.random.choice(N-1, b, replace=False) to draw $b$ different indices of contacts, loop over each infectious individual $I$. This method becomes slow when the infectious population size gets large. After thinking on this, we found we can consider from the susceptible individual's perspective.

\begin{aligned} p_1 &= \text{Pr}\ (\text{an }S \text{ is not contacted by an }I) \\ &= \frac{N-1-b}{N-1} \\ p_2 &= \text{Pr}\ (\text{an }S \text{ is not contacted by all }I \text{'s}) \\ &= p_1^{\vert I \vert} \\ p_3 &= \text{Pr}\ (\text{an }S \text{ is contacted by any }I) \\ &= 1-p_2 \end{aligned}

• Hence we can use np.random.choice to draw $p_3 \cdot \vert S\vert$ individuals from the index set of $S$, without any loops. The speed is 970x faster when $N=10000$. The computation times for the above three methods are shown in the log-log plot below

• Finally, in the spatial case, we need to call a tree query for each $I$ to find its susceptible neighbors, which involes $\vert I \vert$ iterations. This is inevitable since the interaction is determined by the distance. Fortunately, we can again consider from the susceptible individual's perspective, to check if there is any $I$ in an $S$'s neighborhood. More specifically, in the early period, $\vert I \vert \ll \vert S \vert$, so query by $I$ is more efficient. After a time point, $\vert I \vert > \vert S \vert$, we can then change to query by $S$. The complexity is reduced from $O(\vert I \vert)$ to $O(\min(\vert I \vert,\vert S \vert)$. This effectively facilitates our simulation.

PDE model:

• we break the 3 PDEs into 3 systems of $M^2$ ODEs, so when $M$ is large, solving the ODE can be really slow.

• computational performance also heavily depends on $p$. When $p$ is small (of the order $10^{-4}$ or smaller), the ODEs can be solved in a matter of seconds, but when $p$ becomes large (especially larger than $0.01$), the ODEs can take minutes to solve. We haven't completely understand why $p$ has such a big effect on computation, because the interactions of the ODEs are complicated and neither of us has experiences on PDE.

3.3 Limitation

Agent-based model:

• we assumed random direction of movement. In reality, the direction of movement, especially when the movement is relative long distance, is often limited by travel availabilities such as existence of trains, planes, etc.

• we assumed a single rate of movement, $p$. In reality, the rate of movement often changes depending on the destination. For example, if the disease breaks out at the corner, then most rational people would want to avoid going to that corner, thus reducing their rate of movement towards that direction.

• we assumed uniform population density. In reality, urban areas and rural areas are very different in terms of population densities and thus the potential to spread the disease.

PDE model:

• we used second order finite differencing to approximate the Laplacian and changed the PDE into a system of ODEs. This method can generate errors, especially on the corners and the sides. This methods also requires a well-chosen step size, or $M$, because when $M$ is too large, the system of ODEs becomes intractable to solve as we break each PDE into $M^2$ ODEs, but when $M$ is too small, the step size becomes too large and thus generates large errors.

3.4 Further Investigation

We can find directions of further investigation from the limitations discussed above.

Agent-based model:

• we can model the population with the unit square as a graph, where the nodes are the urban areas and the edges are the available traveling routes connecting the urban areas, such as highways, railroads, etc.

• we can limit the movement of individuals to one of the neighboring towns via one of the traveling routes.

• we can change the rate of movement $p$ depending on the travel availability and the destination. For example, $p$ would be large if there are planes connecting two cities and would be small if there's only highways; $p$ would be small, if not zero, if the destination has high rate of infection and would be large if the destination is relatively safe.

PDE model:

• we can try to solve the PDE directly, either using theoretical derivations, or numerical methods, without changing it into a huge system of ODEs.

• we can give $s$, $i$, and $r$ different diffusion parameters instead of using a single parameter $p$. This is because on top of $b$ and $k$, we would still expect some different rate of diffusion for the three populations. For example, if vaccines are invented, people would be taking them at a much faster rate than the disease spread in the first place, so $r$ would have a large diffusion parameter than $i$. Of course, ideally these diffusion parameters, just like $b$ and $k$ as discussed above in the first variation, should be also be able to change with time.

4. Bibliography

1. Project Requirements
2. SEIR Model
3. US Covid-19 Data from CDC