Search

In the infectious disease literature, significant effort has been devoted to studying dynamics at a single scale. For example, compartmental models describing population-level dynamics are often formulated using differential equations. In cases where small numbers or noise play a crucial role, these differential equations are replaced with memoryless Markovian models, where discrete individuals can be members of a compartment and transition stochastically. Classic stochastic simulation algorithms, such as Gillespie's algorithm and the next reaction method, can be employed to solve these Markovian models exactly. The intricate coupling between models at different scales underscores the importance of multiscale modelling in infectious diseases. However, several computational challenges arise when the multiscale model becomes non-Markovian. In this paper, we address these challenges by developing a novel exact stochastic simulation algorithm. We apply it to a showcase multiscale system where all individuals share the same deterministic within-host model while the population-level dynamics are governed by a stochastic formulation. We demonstrate that as long as the within-host information is harvested at a reasonable resolution, the novel algorithm we develop will always be accurate. Moreover, the novel algorithm we develop is general and can be easily applied to other multiscale models in (or outside) the realm of infectious diseases., Comment: 23 pages, 7 figures

Many multi-scale systems exhibit stochastic and deterministic dynamics. Capturing these aspects in a compartmental model is challenging. Notably, low occupancy compartments exhibit stochastic dynamics and high occupancy compartments exhibit deterministic dynamics. Failing to account for stochasticity in small populations can produce 'atto-foxes', e.g. in the classic Lotka-Volterra ODE system. This can be resolved with discrete stochastic models e.g. continuous time Markov chains (CTMCs). Unfortunately, simulating CTMCs is infeasible for most realistic populations. We develop a hybrid CTMC-ODE approach, 'Jump-Switch-Flow' (JSF), available as a Python package. JSF couples CTMCs and ODEs, and depending on the state of the system, individual compartments 'switch' between stochastic 'jumping' and deterministic 'flowing' representations. Consequently, compartments can reach extinct states, thereby resolving atto-type problems. JSF has the desired behaviours of exact CTMC simulation, but is substantially faster that existing alternatives. JSF's value for inference problems, particularly multi-scale problems, is demonstrated with a case-study using a particle filter. With existing data of SARS-CoV-2 infections, we perform a novel analysis to quantify within-host viral clearance. JSF offers a novel approach to compartmental model development and simulation.

The epidemiological behavior of Plasmodium vivax malaria occurs across spatial scales including within-host, population, and metapopulation levels. On the within-host scale, P. vivax sporozoites inoculated in a host may form latent hypnozoites, the activation of which drives secondary infections and accounts for a large proportion of P. vivax illness; on the metapopulation level, the coupled human-vector dynamics characteristic of the population level are further complicated by the migration of human populations across patches with different malaria forces of (re-)infection. To explore the interplay of all three scales in a single two-patch model of Plasmodium vivax dynamics, we construct and study a system of eight integro-differential equations with periodic forcing (arising from the single-frequency sinusoidal movement of a human sub-population). Under the numerically-informed ansatz that the limiting solutions to the system are closely bounded by sinusoidal ones for certain regions of parameter space, we derive a single nonlinear equation from which all approximate limiting solutions may be drawn, and devise necessary and sufficient conditions for the equation to have only a disease-free solution. Our results illustrate the impact of movement on P. vivax transmission and suggest a need to focus vector control efforts on forest mosquito populations. The three-scale model introduced here provides a more comprehensive framework for studying the clinical, behavioral, and geographical factors underlying P. vivax malaria endemicity., Comment: 42 pages and 9 figures

Enzyme-catalysed reactions involve two distinct timescales. There is a short timescale on which enzymes bind to substrate molecules to produce bound complexes, and a comparatively long timescale on which the complex is transformed into a product. The rate at which the substrate is converted into product is characteristically non-linear and is traditionally derived by applying singular perturbation theory to the system's governing equations. Central to this analysis is the assumption that complex formation is effectively instantaneous on the timescale over which significant substrate degradation occurs. This prevents accurate modelling of enzyme kinetics by many particle-based simulations of reaction-diffusion systems as they rely on proximity-based reaction conditions that do not correctly model the fast reactions associated with the complex on the long timescale. In this paper we derive a new proximity-based reaction condition that correctly incorporates the reactions that occur on the short timescale for a specific enzymatic system. We present proof of concept particle-based simulations and demonstrate that non-linear reaction rates typical of enzyme kinetics can be reproduced without needing to explicitly simulate reactions on the short timescale., Comment: 18 pages, 6 figures

\textit{Plasmodium vivax} is the most geographically widespread malaria-causing parasite resulting in significant associated global morbidity and mortality. One of the factors driving this widespread phenomenon is the ability of the parasites to remain dormant in the liver. Known as hypnozoites, they reside in the liver following an initial exposure, before activating later to cause further infections, referred to as relapses. As around 79-96$\%$ of infections are attributed to relapses, we expect it will be highly impactful to apply treatment to target the hypnozoite reservoir to eliminate \textit{P. vivax}. Treatment with a radical cure to target the hypnozoite reservoir is a potential tool to control or eliminate \textit{P. vivax}. We have developed a multiscale mathematical model as a system of integro-differential equations that captures the complex dynamics of \textit{P. vivax} hypnozoites and the effect of hypnozoite relapse on disease transmission. Here, we use our model to study the anticipated effect of radical cure treatment administered via a mass drug administration (MDA) program. We implement multiple rounds of MDA with a fixed interval between rounds, starting from different steady-state disease prevalences. We then construct an optimisation model to obtain the optimal MDA interval. We also incorporate mosquito seasonality in our model to study its effect on the optimal treatment regime. We find that the effect of MDA interventions is temporary and depends on the pre-intervention disease prevalence (and choice of model parameters) as well as the number of MDA rounds under consideration. We find radical cure alone may not be enough to lead to \textit{P. vivax} elimination under our mathematical model (and choice of model parameters) since the prevalence of infection eventually returns to pre-MDA levels., Comment: 39 pages, 14 figures

Turing patterns in reaction-diffusion (RD) systems have classically been studied only in RD systems which do not explicitly depend on independent variables such as space. In practise, many systems for which Turing patterning is important are not homogeneous with ideal boundary conditions. In heterogeneous systems with stable steady states, the steady states are also necessarily heterogeneous which is problematic for applying the classical analysis. Whilst there has been some work done to extend Turing analysis to some heterogeneous systems, for many systems it is still difficult to determine if a stable patterned state is driven purely by system heterogeneity or if a Turing instability is playing a role. In this work, we try to define a framework which uses numerical continuation to map heterogeneous RD systems onto a sensible nearby homogeneous system. This framework may be used for discussing the role of Turing instabilities in establishing patterns in heterogeneous RD systems. We study the Schnakenberg and Gierer-Meinhardt models with spatially heterogeneous production as test problems. It is shown that for sufficiently large system heterogeneity (large amplitude spatial variations in morphogen production) it is possible that Turing-patterned and base states become coincident and therefore impossible to distinguish. Other exotic behaviour is also shown to be possible. We also study a novel scenario in which morphogen is produced locally at levels that could support Turing patterning but on intervals/patches which are on the scale of classical critical domain lengths. Without classical domain boundaries, Turing patterns are allowed to bleed through; an effect noted by other authors. In this case, this phenomena effectively changes the critical domain length. Indeed, we even note that this phenomena may also effectively couple local patches together and drive instability in this way., Comment: 10 figures

Malaria is an infectious disease with an immense global health burden. Plasmodium vivax is the most geographically widespread species of malaria. Relapsing infections, caused by the activation of liver-stage parasites known as hypnozoites, are a critical feature of the epidemiology of Plasmodium vivax. Hypnozoites remain dormant in the liver for weeks or months after inoculation, but cause relapsing infections upon activation. Here, we introduce a dynamic probability model of the activation-clearance process governing both potential relapses and the size of the hypnozoite reservoir. We begin by modelling activation-clearance dynamics for a single hypnozoite using a continuous-time Markov chain. We then extend our analysis to consider activation-clearance dynamics for a single mosquito bite, which can simultaneously establish multiple hypnozoites, under the assumption of independent hypnozoite behaviour. We derive analytic expressions for the time to first relapse and the time to hypnozoite clearance for mosquito bites establishing variable numbers of hypnozoites, both of which are quantities of epidemiological significance. Our results extend those in the literature, which were limited due to an assumption of non-independence. Our within-host model can be embedded readily in multi-scale models and epidemiological frameworks, with analytic solutions increasing the tractability of statistical inference and analysis. Our work therefore provides a foundation for further work on immune development and epidemiological-scale analysis, both of which are important for achieving the goal of malaria elimination., Comment: This is a pre-copyedit version of an article published in Bulletin of Mathematical Biology. The final authenticated version is available online at: https://doi.org/10.1007/s11538-020-00706-1

In 1917, Marian von Smoluchowski presented a simple mathematical description of diffusion-controlled reactions on the scale of individual molecules. His model postulated that a reaction would occur when two reactants were sufficiently close and, more specifically, presented a succinct relationship between the relative proximity of two reactants at the moment of reaction and the macroscopic reaction rate. Over the last century, Smoluchowski reaction theory has been applied widely in the physical, chemical, environmental and, more recently, the biological sciences. Despite the widespread utility of the Smoluchowski theory, it only describes the rates of second order reactions and is inadequate for the description of higher order reactions for which there is no equivalent method for theoretical investigation. In this paper, we derive a generalised Smoluchowski framework in which we define what should be meant by proximity in this context when more than two reactants are involved. We derive the relationship between the macroscopic reaction rate and the critical proximity at which a reaction occurs for higher order reactions. Using this theoretical framework and using numerical experiments we explore various peculiar properties of multimolecular diffusion-controlled reactions which, due to there being no other numerical method of this nature, have not been previous reported.

Spatial reaction-diffusion models have been employed to describe many emergent phenomena in biological systems. The modelling technique most commonly adopted in the literature implements systems of partial differential equations (PDEs), which assumes there are sufficient densities of particles that a continuum approximation is valid. However, due to recent advances in computational power, the simulation, and therefore postulation, of computationally intensive individual-based models has become a popular way to investigate the effects of noise in reaction-diffusion systems in which regions of low copy numbers exist. The stochastic models with which we shall be concerned in this manuscript are referred to as `compartment-based'. These models are characterised by a discretisation of the computational domain into a grid/lattice of `compartments'. Within each compartment particles are assumed to be well-mixed and are permitted to react with other particles within their compartment or to transfer between neighbouring compartments. We develop two hybrid algorithms in which a PDE is coupled to a compartment-based model. Rather than attempting to balance average fluxes, our algorithms answer a more fundamental question: `how are individual particles transported between the vastly different model descriptions?' First, we present an algorithm derived by carefully re-defining the continuous PDE concentration as a probability distribution. Whilst this first algorithm shows strong convergence to analytic solutions of test problems, it can be cumbersome to simulate. Our second algorithm is a simplified and more efficient implementation of the first, it is derived in the continuum limit over the PDE region alone. We test our hybrid methods for functionality and accuracy in a variety of different scenarios by comparing the averaged simulations to analytic solutions of PDEs for mean concentrations., Comment: MAIN - 24 pages, 10 figures, 1 supplementary file - 3 pages, 2 figures

In this paper, three multiscale methods for coupling of mesoscopic (compartment-based) and microscopic (molecular-based) stochastic reaction-diffusion simulations are investigated. Two of the three methods that will be discussed in detail have been previously reported in the literature; the two-regime method (TRM) and the compartment-placement method (CPM). The third method that is introduced and analysed in this paper is the ghost cell method (GCM). Presented is a comparison of sources of error. The convergent properties of this error are studied as the time step $\Delta t$ (for updating the molecular-based part of the model) approaches zero. It is found that the error behaviour depends on another fundamental computational parameter $h$, the compartment size in the mesoscopic part of the model. Two important limiting cases, which appear in applications, are considered: (i) \Delta t approaches 0 and h is fixed; and (ii) \Delta t approaches 0 and h approaches 0 such that \Delta t/h^2 is fixed. The error for previously developed approaches (the TRM and CPM) converges to zero only in the limiting case (ii), but not in case (i). It is shown that the error of the GCM converges in the limiting case (i). Thus the GCM is superior to previous coupling techniques if the mesoscopic description is much coarser than the microscopic part of the model.

The two-regime method (TRM) has been recently developed for optimizing stochastic reaction-diffusion simulations. It is a multiscale (hybrid) algorithm which uses stochastic reaction-diffusion models with different levels of detail in different parts of the computational domain. The coupling condition on the interface between different modelling regimes of the TRM was previously derived for one-dimensional models. In this paper, the TRM is generalized to higher dimensional reaction-diffusion systems. Coupling Brownian dynamics models with compartment-based models on regular (square) two-dimensional lattices is studied in detail. In this case, the interface between different modelling regimes contain either flat parts or right-angled corners. Both cases are studied in the paper. For flat interfaces, it is shown that the one-dimensional theory can be used along the line perpendicular to the TRM interface. In the direction tangential to the interface, two choices of the TRM parameters are presented. Their applicability depends on the compartment size and the time step used in the molecular-based regime. The two-dimensional generalization of the TRM is also discussed in the case of corners.

Two algorithms that combine Brownian dynamics (BD) simulations with mean-field partial differential equations (PDEs) are presented. This PDE-assisted Brownian dynamics (PBD) methodology provides exact particle tracking data in parts of the domain, whilst making use of a mean-field reaction-diffusion PDE description elsewhere. The first PBD algorithm couples BD simulations with PDEs by randomly creating new particles close to the interface which partitions the domain and by reincorporating particles into the continuum PDE-description when they cross the interface. The second PBD algorithm introduces an overlap region, where both descriptions exist in parallel. It is shown that to accurately compute variances using the PBD simulation requires the overlap region. Advantages of both PBD approaches are discussed and illustrative numerical examples are presented., Comment: submitted to SIAM Journal on Applied Mathematics