Applying the chemical-reaction definition of mass action to infectious disease modelling

The law of mass action is used to govern interactions between susceptible and infected individuals in a variety of infectious disease models. However, the commonly used version is a simplification of the version originally used to describe chemical reactions. We reformulate a general disease model using the chemical-reaction definition of mass action incorporating both an altered transmission term and an altered recovery term in the form of positive exponents. We examine the long-term outcome as these exponents vary. For many scenarios, the reproduction number is either 0 or $\infty$, while it obtains finite values only for certain combinations. We found conditions under which endemic equilibria exist and are unique for a variety of possible exponents. We also determined circumstances under which backward bifurcations are possible or do not occur. The simplified form of mass action may be masking generalised behaviour that may result in practice if these exponents ``fluctuate'' around 1. This may lead to a loss of predictability in some models.


Introduction
The law of mass action is a standard modelling tool that has been used to characterize the transmission of disease between infected and susceptible individuals. As such, the law of mass action has had far-reaching implications in the prediction of disease trajectories and the evaluation of health interventions.
Initially, however, the law of mass action was first applied to describe the reaction rate between chemical reactants [20], portraying the process by which chemicals collide and interact to form different chemical combinations. The law of mass action under the context of chemical reactions is well known [19] and has been generalized to include time-dependent rate constants, real-valued kinetic orders [9,17,7] and even power-law kinetics to describe multi-chemical reactions [14,15,16]. Importantly, the use of power-law kinetics to describe multi-chemical reactions permits the possibility of the complete exhaustion of chemical reactants, thereby potentially achieving a steady state solution in finite time [16].
The law of mass action in the infectious-disease context is also well studied. Traditionally in infectious-disease modelling, the law of mass action is often simplified from the original formulation proposed by Guldberg and Waage [20] to a bilinear function of susceptible and infected proportions multiplied by a transmission rate. Generalizations of the law of mass action do exist, however, and are usually developed to account for behavioural characteristics, like the crowding of infected individuals, measures taken by susceptible individuals to avoid exposure to infection or the effects of behavioural changes [10,11,18,13,1]. Generalizations of the law of mass action to power-law kinetics have also been considered [10,11]; however, due to the assumption of linearity often applied to the recovery term, the complete exhaustion (or elimination) of the infected proportion is typically not physically possible. That said, there does exist substantial motivation for using a nonlinear recovery rate that may permit the exhaustion (or elimination) of the infected proportion. For instance, a nonlinear recovery rate could account for additional behavioural characteristics, like active treatment that includes quarantine or isolation [22,21], the effect of limited treatment resources [3,12] or the effects of population demographics [6].
There have been a handful of infectious-disease models that have used variations of the generalised mass action transmission function. Wang et al. used an infection form with α = 2 to represent the effect of double exposure of susceptibles in a model with medical-resource contraints [4]. Ermentrout & Ermentrout [5] used the same infection form to represent two zombies needed to infect one human. Beeton et al. [2] used a more generalised form (with powers in the infection term but also the susceptible term) to model the effects of stochasticity in a zombie model. Greenhalgh et al. used this form to describe measles in Iceland and its re-emergence from Denmark via fishermen [6]. The use of exponents in models has been used to quantify activities of cell production and removal in bone remodelling [8].
Here we quantify the effects of generalizing both the law of mass action and the recovery rate in infectious-disease-transmission dynamics. In particular, we examine the effects of power-law kinetics, as they can be used to generalize both the law of mass action and the recovery rate. Our results include power-law-exponent conditions for stability, the existence of endemic equilibria and identifying backward bifurcations.
We address the following research questions: Can we quantify the effect of the exponents in the model? Specifically, can we calculate R 0 and determine stability for various choices of the exponents? Can we determine the existence of endemic equilibria or conditions where they do not arise? Does the generalised model lead to backward bifurcations?

The Model
The general form of the mass-action transmission model is given by where S, I and R are proportions of susceptible, infected and recovered individuals. The population is assumed to be growth neutral so that the birth/death rates are given by b. The infection rate is β and the recovery rate is η. Note that all parameters are assumed to be positive. The exponents α and ξ describe a modified form of infection and recovery. Here, we generalise α to take values from the positive real numbers. The exponent ξ may account for differences in recovery related to the immune system. Different people may recover faster or slower, depending on how effective their immune systems are. The usual form of mass-action transmission is the case when α = ξ = 1. We can assume α, ξ > 0. The equations (2.1)-(2.3) are supplemented with the initial conditions S(0) = S 0 , I(0) = I 0 , R(0) = R 0 . (2.4) We will also assume that the system (2.1)-(2.4) is epidemiologically well defined; i.e., (S 0 , I 0 , R 0 ) ∈ [0, 1] 3 such that S 0 = 0 and I 0 = 0 (2.5) Note that the model may not have a unique solution when α or ξ is less than one, as the right-hand side of equation (2.2) will not be Lipschitz around I = 0. Also, there is a possibility that no global solution exists for the model when α or ξ is greater than one. Proof.
(1) Existence of solutions. We define It is clear that the right-hand side of the system (2.1)- With S 0 and I 0 satisfying (2.4) and (2.5), we have S(t), I(t) > 0 for t ≥ 0. Starting with which implies Integrating on [0, t) , t > 0, gives N (t) = 1 for all t > 0. Since X(0) ≥ 0 for all t > 0, we deduce that X(t) ≤ 1 and thus bounded.
Remark. Theorem 3.1 proves local solutions exist, although not necessarily global.

3.2.
Stability of the disease-free equilibrium. The disease-free equilibrium is (S, I, R) = (1, 0, 0). To examine stability, we set S = 1 and then have Setting the part in the parentheses equal to zero and rearranging, we have We have several cases: This is the value for the standard model.
This is summarised in Figure 1. An important observation is how sensitive the standard value of R 0 is to fluctuations in α and ξ when these values are close to 1. Since the disease will invade if R α,ξ 0 > 1 and will be eradicated if R α,ξ 0 < 1, it follows that small changes in either exponent can have drastic effects on the outcome. This is particularly true when the exponents are both close to 1 (as in the standard model).
Substitution into the S = 0 equation yields It follows that the existence of any endemic equilibria (S * , , depends on α and ξ.
We set In the following proposition, we provide sufficient conditions for the existence of the EE.
Proof. We introduce the function which is an antiderivative of f (I) give in (3.3). Using the mean-value theorem for 0 ≤ I ≤ 1, there exists at least one 0 (3.6) Note that both f (0), f (1) < 0 and f has a unique local maximum at the critical value Figure 2). (b) ξ = 1. Assume that β ≥ β c > 0. Then, using Proposition 3.2, there exists at least one EE.
No conclusion can be drawn from applying the IVT on f . If β ≥ β c > 0, then Proposition 3.2 implies the existence of at least one EE. This case is shown in Figure 10.
1) < 0 and we have no conclusion from f . Assuming that β ≥ β c > 0, then, using Proposition 3.2, there exists at least one EE. This case is shown in Figure 12.

Note that
Now, if R α,ξ 0 = ∞, then the DFE is unstable. We thus need to consider only the following cases.
One unique EE for β > η No conclusion for β < η Possible At least one EE for R α,ξ From the EE, we have For ξ = 2 we have Hence the unique EE is globally stable for α = 1 and ξ = 2, when R 1,2 0 > 1.

Discussion
We modified an SIR mass-action model to include I α and I ξ with α, ξ > 0 in the infection and recovery terms, respectively. The powers α and ξ account for differences in infectivity and recovery. We used α and ξ to calculate the basic reproduction number. Small changes in either of the exponents change R 0 between 0, different finite values and ∞. As a result, the stability of the DFE is sensitive to the choice of α and ξ and requires further restrictions to have a globally stable DFE. Similar results show that the existence of the EE is also sensitive to the choice of the exponents. Numerical simulations are provided for the time-series solution for different values of α and ξ to illustrate our theoretical results.
Since the disease will invade if R α,ξ 0 > 1 and can be eradicated if R α,ξ 0 < 1, it follows that small changes in either exponent can have drastic effects on the outcome. This is particularly true when the exponents are both close to 1 (as in the standard model). It follows that epidemic models that use the original mass-action formulation may not have meaningful predictions if a reproduction number is calculated.
Some limitations apply, which should be acknowledged. We assumed that the population growth was neutral; that is, rates of growth and death are equal (b), which is not true in general. Mass-action transmission, whether in its original form or the form used in most epidemic models, assumes that populations are well mixed, an assumption based on the random mixing of particles from which the original term was derived. Such an assumption rarely applies in practice.
Future work will apply these results to standard incidence. While we were guided by the chemicalreaction definition of mass action, it would be interesting to generalise all other terms in the model. The rich behaviour exhibited here suggests that other extensions would like yield a range of unanticipated outcomes. Likewise, although different forms of the transmission function have been investigated, few papers have looked at different forms of the recovery rate. Extending this beyond our formulation would doubtless result in further unexpected phenomena.
Given the sensitivity of our results to small changes in α and ξ, especially close to α = ξ = 1 (the standard case), it is likely that many disease models are sitting on a knife-edge situation. If such fluctuations occur in the real world, then the ability of disease models to make long-term predictions may be vastly curtailed.