On the use of reactive multiparticle collision dynamics to gather particulate level information from simulations of epidemic models

Transmission of a communicable disease is sensitive to person-to-person contacts. The contact patterns can help identify people at risk of contracting the disease and where an outbreak can be prevented. Modern epidemiology relies heavily on mathematical modeling to comprehend the transmission dynamics of an infectious disease. The knowledge gained from mathematical models is vital for assessing the impact as well as improving the effectiveness of intervention strategies such as lockdowns, vaccination, social distancing, etc., in controlling the spread.

The most common formalism for the mathematical description of disease transmission is the deterministic or compartmental approach, which assumes that the sizes of the compartments are large enough that the mixing of individuals is homogeneous, and thus, the law of mass action holds. However, at the onset of a disease outbreak, most individuals in a population are susceptible and the number of infectives is small. The infection transmission depends largely on the contact pattern between the individuals. The mass action assumption should therefore be replaced by a model that takes into account stochastic effects.¹ 

From a kinetic theory standpoint, an epidemic model can be viewed as a chemically reacting system of species that follows a continuous-time Markov chain dynamics with discrete state space.² Individuals are grouped according to their state, with commonly used states being “susceptible (S),” “infectious (I),” and “recovered (R).” These states correspond to the chemical species in such an analog model. The state of the system at time t is the random vector containing the discrete number of individuals, or molecules, of each species. Reaction rates are then used to incorporate all transmission dynamics into the epidemic model.

In recent years, there has been a growing interest in developing methods of statistical mechanics to model a population as a set of interacting particles with spatial dynamics. These methods have also been applied to describe the evolution and spread of epidemics. One such contribution was by Delitala,³ who proposed epidemic modeling in the framework of kinetic (Boltzmann) models. The models take into account stochastic type interactions between individuals and describe the space–time evolution of healthy and infected individuals. The application of Boltzmann-type kinetic models toward epidemiology can also be found in Refs. 4–8. The models use thermodynamic properties (position, velocity, temperature, etc.) to describe the interaction dynamics in a similar way to elastic collisions between fluid particles, naturally introducing noise/randomness in the system. It is worth noting that there are situations, such as people moving in a shopping mall or a supermarket, and children moving in a classroom or school campus, where the movement is altered by ambient noise, and thus, it is plausible to regard individuals as random particles.^8,9

A kinetic Monte Carlo simulation algorithm was presented in Ref. 10 to test different lockdown patterns to monitor the spread of novel coronavirus. The method is based on modeling a population on the principles of Brownian motion, and the population characteristics are described using thermodynamic parameters (temperature, density, and diffusion). The time period when there is no lockdown is characterized by high temperature, whereas a lockdown is characterized by low temperature. Results show that social distancing and avoiding creation of hotspots are effective strategies to control the infection rate over a long period.

The application of Brownian Dynamics (BD) to a susceptible-infected-recovered (SIR) model was discussed in Ref. 9. Effects of contagion rate, particle density, and particle velocity are discussed. The results are in good agreement with the traditional model. Critical density (i.e., the density above which the individuals are immediately infected) and contagious radius, both of which facilitate the spread of virus, are explored. In addition, an analytical expression for the contagion rate in terms of microscopic parameters is obtained using a mean-free path analysis, which is a widely used concept in the kinetic theory of gases.

In Ref. 11, a susceptible-exposed-infectious-recovered (SEIR) model was analyzed using the molecular dynamics simulation method. Although molecular dynamics is a useful technique in statistical mechanics, it is extremely time consuming. An efficient alternative is the Direct Simulation Monte Carlo (DSMC) method,¹² which is typically used to study rarefied gas flow problems^13,14 and has also found interesting applications in physics and chemistry.^15,16 Recently, Pacheco et al.¹⁷ investigated the use of the DSMC method to model the COVID-19 pandemic. The model is validated by fitting it to COVID-19 case counts in New York City. Sensitivity analysis showed that the model is sensitive to the number of initially infected individuals, interaction diameter (i.e., the region surrounding an infected individual within which a passerby has a high probability of catching the infection), temperature (i.e., mobility of individuals), and disease progression rate. The model suggests that while positive cases, hospitalizations, and deaths are important indicators for “re-opening,” understanding how the asymptomatic population is trending is vital. A high number of asymptomatic cases can result in spikes even if fewer hospitalizations or deaths are recorded.

One of the cornerstones of statistical mechanics is dynamical density functional theory (DDFT), which is an extension of classical density functional theory (DFT) to nonequilibrium systems. Originally developed for the study of fluids, DDFT now has diverse applications. A comprehensive review of DDFT can be found in Ref. 18. In Ref. 19, DDFT was used to model the effects of social distancing and isolation in preventing the spread of COVID-19. A DDFT-based extension of an SIR model is presented. The model considers individuals as diffusing particles with repulsive interactions that correspond to social distancing and isolation. These repulsive interactions are shown to reduce the number of infections as well as allow flattening of the curve.

De Rosis²⁰ simulated epidemic dynamics using the lattice Boltzmann method (LBM) and demonstrated excellent agreement between numerical results and the solution of the SIR model with spatial effects. Two advantages of the LBM are worth mentioning. First, the reactive part is introduced by simply adding a source term projected through the weights into the distribution space. Second, unlike conventional SIR models with diffusion, computation of spatial derivatives is not required in the LBM formulation since spatial effects are naturally accounted for in the nonreactive part, i.e., the collisions. A limitation of the LBM is that it tends to be unstable for reduced diffusivity and an increase in contact and recovery rates, as illustrated by linear stability analysis.

A study of spatiotemporal dynamics of a particle-based SIR model was presented in Ref. 21. The model is simulated using a lattice gas cellular automaton (LGCA), a precursor of the LBM.²² Effects of spatial inhomogeneities on the epidemic dynamics and vaccination strategies are investigated. In addition, approximate mean-field equations for the automaton are derived. Comparison of simulation results from the mean-field description with those obtained from the LGCA approach demonstrates that the latter is a potential alternative to the macroscopic SIR model with diffusion. Other important studies related to particle-based modeling for spread of infectious disease include Refs. 23–30.

A recent particle-based method called reactive multiparticle collision (RMPC) dynamics^31,32 is currently a popular technique for simulating spatially distributed chemically reacting systems.^33–36 The method has not been applied to epidemic models, and its potential in this area remains unexplored. The aim of this paper is to investigate the performance of RMPC in simulating an epidemic model and obtaining useful epidemiological statistics such as average contact duration, average initial infection time, etc., which cannot be determined via traditional compartmental models or population-based stochastic models.

The paper is outlined as follows: Sec. II describes the details of RMPC dynamics. In Sec. III, we present results on contact information, reinfection statistics, and the timing of exposure obtained from RMPC simulations of an SIR, SIRS, and SEIR model. Conclusions to this paper are given in Sec. IV.

Reactive multiparticle collision (RMPC) dynamics is the extension of multiparticle collision (MPC) dynamics³⁷ for diffusive motion to include reactions. RMPC describes the particle dynamics by three steps: a streaming step, a collision step, and a reaction step. Consider a system of N point particles in a volume V, which is divided into N_c cells of volume V_c labeled by indices ξ. Each cell is assigned a random rotation operator $ω̂ξ$ chosen from a set Ω. Suppose the ith particle has position r_i and velocity v_i at time t. The particles undergo local collisions and reactions at discrete time steps τ followed by free-streaming.

Disease transmission is often facilitated by social contacts, and therefore, contact information is essential to controlling the transmission. Here, we demonstrate that RMPC is able to determine the number of contacts and the duration of the contacts by simulating an SIR model.

The model is simulated in the computational domain 1 × 1 × 1 km³ that is divided into 20 × 20 × 20 cubic cells of side length h = 0.05 km subject to periodic boundary conditions and k_BT/m = 75 km²/day². We choose N = 40, 000 with S(0) = 39, 800, I(0) = 200 and R(0) = 0. The parameters β and γ are set as 1.0 and 0.2 day⁻¹, respectively. We run a single realization of RMPC dynamics. To obtain contact data with high temporal resolution, RMPC simulation is performed on a second-by-second basis (that is, τ = 1 s) for a day. This results in a mean-free path, $λ=kBT/mτ≈0.0001$ km, smaller than the cell size h. We define a contact as a susceptible being present in the same collision cell (e.g., room, office, etc.) as an infectious individual. This is relevant for the spread of airborne infections. A distribution of the number of contacts is given in Fig. 3(a). We also determined the duration of contacts, i.e., the total time a susceptible remained in the same cell as an infectious person, and Fig. 3(b) shows the distribution of the contact duration.

Our simulation shows that number of contacts and the average contact duration (total duration of contacts divided by the total number of contacts) are sensitive to changes in the transmission rate β and average population speed. That is to say, when the transmission rate β increases, the number of contacts increases [Fig. 4(a)] whereas the average contact duration decreases [Fig. 4(b)]. Similarly, the number of contacts increases and the average contact duration decreases with increasing values of average population speed, as shown in Figs. 5(a) and 5(b), respectively. This observation indicates that at lower speed, a longer contact duration will result in an epidemic spread. On the other hand, when individual mobility is high, infections can occur from frequent short-duration contacts.

It is assumed in the SIR model that a person is completely immune upon recovery from the disease. When reinfections are possible, as in the case of COVID-19, an SIRS model is considered. Recovered individuals become susceptible at a rate ɛ, as shown in Fig. 6. The transition rates for the model are shown in Table II.

For parameters, we use β = 3/14 day⁻¹, γ = 1/14 day⁻¹, and ɛ = 1/180 day⁻¹. We simulate this model using RMPC, with reflective boundary conditions and k_BT/m = 2 km²/day², in the 1 km × 1 km × 0.2 km domain that is represented as 5 × 5 cells of side length h = 0.2 km. Here, N = 200 with S(0) = 199, I(0) = 1, and R(0) = 0. A single realization of RMPC is performed over a period of 1 year using τ = 0.1 day, leading to the mean-free path λ = 0.141 42 km, which is slightly smaller than the cell size h. We obtained the temporal evolution of S, I, and R populations [Fig. 7(a)] and the spatial dynamics of the epidemic [Figs. 7(b)–7(f)]. The latter also show an individual being tracked who was initially susceptible, caught infection for the first time on Day 8, recovered on Day 16, became susceptible again on Day 63, and was reinfected on Day 91. For all such individuals, we recorded the first infection time [Fig. 8(a)] and the time of reinfection [Fig. 8(b)]. We found that the average time that an individual first becomes infected (30.77 days) is much faster than becoming reinfected (98.60 days).

Most infectious diseases have a latent period, which is the time elapsed between being infected to being infectious or contagious. The simplest way to incorporate disease latency in a mathematical model is by grouping such individuals in compartment E (for exposed), as shown in Fig. 9. Upon being infected, the individuals move to this compartment at a rate βSI/N and remain there for a period of 1/σ (the mean latent period) before moving to the I compartment at a rate σE (Table III).

We choose β = 3.0 day⁻¹, σ = 1/8 day⁻¹, and γ = 1/5 day⁻¹. The spatial domain and the boundary conditions are the same as for the previous model. Here, N = S + E + I + R with S(0) = 199 and I(0) = 1. Using τ = 0.1 day, we obtained temporal dynamics of the four sub-populations [Fig. 10(a)] and dynamic spreading of infection on Day 7 and Day 14 [Figs. 10(b) and 10(c)] from a single realization of RMPC. For the entire population, we recorded the time a susceptible becomes exposed, an exposed becomes infectious, and an infectious person recovers from the disease [Figs. 10(d)–10(f)]. We found that the average times (in days) for exposure, infection, and recovery are 10.25, 17.80, and 22.57, respectively.

In this paper, we evaluate reactive multiparticle collision (RMPC) dynamics for epidemic modeling. RMPC treats each individual as a particle with position, velocity, and epidemic state. RMPC keeps track of every individual at every time step and is therefore able to determine contact pattern, reinfection statistics, and timing of exposure, infection, and recovery. We demonstrate this by performing RMPC simulations of an SIR, SIRS, and SEIR model to obtain the average contact duration; the average initial infection time and the average re-infection time; and the average times for exposure, infection, and recovery, respectively. Clearly, such detailed/microscopic-level data cannot be obtained from the deterministic solutions of these models or population-based stochastic methods.

RMPC looks promising and can serve as a useful simulation tool when considering more realistic and complex epidemic models. In the future, we would like to study the effectiveness of RMPC by comparing the simulation results with real data. Other possible directions include deducing the speed of propagation of an infection wave in space and quantifying how the effective infection transmission rate is affected by space.

Z. Memon acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) (Funding Reference No. PGS D-570149-2022) and the Ontario Graduate Scholarship. K. Rohlf acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), Discovery Grant No. RGPIN-2017-04672, and 2007 NSERC Research Tools and Infrastructure (RTI) grant. This research has also been partially supported by the Toronto Metropolitan University Mathematics Department and Faculty of Science Dean’s Research Fund. Special thanks go to the anonymous reviewers for their time and input.

The authors have no conflicts to disclose.

Zaib Un Nisa Memon: Conceptualization (lead); Funding acquisition (equal); Investigation (lead); Methodology (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Katrin Rohlf: Conceptualization (supporting); Funding acquisition (equal); Methodology (equal); Software (supporting); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.