Month: March 2020

Physical Modelling of a Pandemic

nafshordi1 Comment

These are strange times, the likes of which we probably see once in a generation, and it is amazing that it is happening across the world simultaneously, and in the age of internet. As physicists, as much as one would like to carry on , it is hard to ignore what is going on around you, and try to understand it, the same way we understand any other physical system. So here it goes:

A Growth+Diffusion Model for Early Epidemic

As we approach any other physics problem, we start from first principles that describe the system of interest. For the purpose of this note, I will start by describing how how a contagion is transmitted within a community and spreads across large geographical locations. What bothers me about the standard treatments (e.g., the SIR equation) is that, by reducing the system to ordinary differential equations, it ignores this geographical diversity and thus (as we shall see below) might miss important phenomenological features of the evolution.

Instead, I will use an inhomogeneous linear growth+diffusion equation to describe the early phase of an epidemic, prior to large-scale immunity or social intervention:

Here n(x,t) is the number density of infected individuals as a function of space and time. R(x) is the local linear growth rate (= ln(2) divided by the doubling time), while D(x) is a diffusion coefficient, quantifying how infected individuals can move around. Since the coefficients don’t have explicit time-dependence, we can decompose the solution into modes with exponential growth (or decay):

where \lambda_m‘s and u_m({\bf x})‘s are eigenvalues and eingenstates of the elliptic operator R+\nabla\cdot D \nabla:

For constant D, this is the same equation as energy states of a particle in 2D that satisfies time-independent Schrodinger Equation:

For example, one result that we can now use is that for generic random potentials (or disorder) in 2D the eigenstates are localized, otherwise known as Anderson localization. One can think of these different localized eigenstates as localized communities with different doubling times for the growth of the epidemic.

The next physics result we shall use is the probability distribution of spacing of energy states of a random potential, known as Wigner surmise:

which we shall use to quantify the distribution of the largest eigenvalue (equivalent to the energy of the most energy bound state). Plugging this into Equation (2) yields:

Now, for large \lambda and t, we would like to use saddle approximation compute this integral. To do this, we have to expand -\ln A(\lambda) to second order in \lambda, i.e. -\ln A(\lambda) = a_0 +a_1 \lambda + a_2 \lambda^2 + \dots where we assume a_2 >-C, otherwise the integral would diverge (or we need to include more terms). The saddle-point approximation for large times yields:

In other words, we predict a super-exponential growth at early times, which is contrast to the exponential expectation from simple uniform models, such as the SIR equation. The reason for this is clearly that we are dealing with a distribution of growth rates, and the later times will be dominated by populations with faster growth rates (smallest doubling time), no matter how small they start.

One clear limitation of Equation (8), apart that it only applies to the linear regime and no interventions, is that at some point we will be dominated by the community with the largest eigenvalue, at which point the exponent should become linear time. This will be followed by further nonlinear effects. I will include these effects, by terms that are higher order in time, i.e.

Comparison to Data

We shall next try to fit this to fatality data for different countries, up until March 29th, 2020 (using data provided by https://ourworldindata.org/coronavirus). We find the best fit parameters:

b_2= (7.7 \pm 1.4) \times 10^{-3}~ {\rm day}^2 and b_3= (8.7 \pm 3.6) \times 10^{-5}~ {\rm day}^3. However, these parameters are highly degenerate, as can be seen in the 1 and 2\sigma contours

Best fit parameters for US (blue), Canada (red) data, and Italy (green) as of March 29, 2020. There is just over 2\sigma evidence for the cubic term in US data, which means that t\exp(b_2 t^2) can describe most of US data. The Canadian parameters are more poorly constrained, but appear perfectly consistent with those of the US. The Italian parameters, are considerably better constrained but still consistent with the two other countries (but see below).

In fact, including the cubic terms only became necessary in the last 3 days to fit the US data. Furthermore, the Canadian parameters are perfectly consistent with those of US, but with larger errors.

Some comments on fitting the data: The most rigorous way to fit the data (assuming that you have a perfect model) is to fit daily values using Poisson statistics, as they are independent. Using this, our best fit to US, Canada, and Italy data gives \chi^2_{\rm red} =1.9, 0.9, and 12 respectively. I interpret this as data not being completely uncorrelated, as expected, since infection happens in clusters. As such, for parameter estimation I divide the Poisson log-likelihood by \chi^2_{\rm red} to reflect this sample variance.

The next figure shows the logarithmic derivative \dot{N}/N for US, Canada, Italy, and South Korea. While for the US and Canada, t=1 is the report of first death, the Italy data is shifted forward by 25 days from the day of first death. The best curve fit is to the US data, but clearly fits both Canada and Italy (after the shift), as the parameters are consistent at 95% confidence level, per the figure above.

I have also added South Korean data, but shifted by 40 days. While the US best-fit, extrapolated forward can fit well the first 20 days of S Korean dat, it clearly misses the second half. However, it could well be that higher order terms in the exponent of Equation (9) will become important at this point.

Logarithmic derivative of the total fatality \frac{d\ln N}{dt} \simeq \frac{N(t)-N(t-1)}{N(t-1)} for US (blue), Canada (red), Italy (green), and South Korea (black), as of March 29th. For US and Canada, t=1 is at the report of first death, while the data for Italy and South Korea are shifted forward by 25 and 40 days, respectively. The curve shows the best fit to US data.

This plot is of course very suggestive. Could there be a universal \dot{N}/N as a function of time (analogous to Hubble constant H(t), or Friedmann equation, in cosmology)? This may suggest the intriguing possibility that the true onset of the Italian (South Korean) epidemic might have been 25 (40) days before the first death was attributed to the Covid-19 epidemic. Of course, the alternative is that there is no universality, and the different behaviors are dictated by environmental and social factors.

Final Note

Finally, let me end on a somber note: The best-fit model (to US, Canada, and Italy) predicts that Canadian mortality will pass 1000 around April 9th, while the US mortality will pass 10,000 around April 5th. I will not dare extrapolate beyond that point, and I sincerely hope that the model is wrong.

So, Stay Home! Stay Safe!