A general model of hormesis in biological systems and its application to pest management

Hormesis, a phenomenon whereby exposure to high levels of stressors is inhibitory but low (mild, sublethal and subtoxic) doses are stimulatory, challenges decision-making in the management of cancer, neurodegenerative diseases, nutrition and ecotoxicology. In the latter, increasing amounts of a pesticide may lead to upsurges rather than declines of pests, ecological paradoxes that are difficult to predict. Using a novel re-formulation of the Ricker population equation, we show how interactions between intervention strengths and dose timings, dose–response functions and intrinsic factors can model such paradoxes and hormesis. A model with three critical parameters revealed hormetic biphasic dose and dose timing responses, either in a J-shape or an inverted U-shape, yielding a homeostatic change or a catastrophic shift and hormetic effects in many parameter regions. Such effects were enhanced by repeated pulses of low-level stimulations within one generation at different dose timings, thereby reducing threshold levels, maximum responses and inhibition. The model provides insights into the complex dynamics of such systems and a methodology for improved experimental design and analysis, with wide-reaching implications for understanding hormetic effects in ecology and in medical and veterinary treatment decision-making. We hypothesized that the dynamics of a discrete generation pest control system can be determined by various three-parameter spaces, some of which reveal the conditions for occurrence of hormesis, and confirmed this by fitting our model to both hormetic data from the literature and to a non-hormetic dataset on pesticidal control of mirid bugs in cotton.


Model formulation
In order to formulate the novel discrete single species model with perturbation within each generation, analytical and piecewise constant methods were used to extend the classic discrete Beverton-Holt and Ricker models [1, 2, 3].

The extended Beverton-Holt model (EBHM)
Consider the following Logistic equation with each integer interval t ∈ [n, n + 1]. Although model (S1.1) can describe the growth of many kinds of populations, including tumor cells and pests, as an example for modelling control tactics we only consider a pest population in the following. Without loss of generality, assuming that there exists a positive constant θ ∈ [0, 1] such that a chemical control tactic was applied at n + θ and a proportion of the pest was killed, denoted by q (characterization of pesticide or drug dose or effectiveness) so that p = 1 − q is the survival rate after the pesticide has been sprayed.
There are two steady states for model (S1.8), which can be calculated as By simple calculation, for the positive equilibrium N * , we have < 0 (S1.9) and dN * dθ = K(Rp − 1)R θ ln(R)(p − 1) [pR − 1 + R θ (1 − p)] 2 < 0. (S1.10) These derivations of N * with respect to p and θ indicate that the later the spraying of pesticide and the larger the killing rate, the more effective is the control in relation to the stable population level, and these results show that the paradox and hormesis could not occur for the Beverton-Holt model when the chemical control action is applied within a discrete generation.

The hormesis Ricker model (HRM)
The piecewise constant method is employed in this subsection to derive the discrete Ricker model with instantaneous perturbations within each generation, i.e., we consider the following piecewise single species model , t ∈ (n, n + θ] (S1.11) with N + θ .
A typical feature of models EBHM and HRM , is that the control measure is applied instantaneously at time point n + θ within two successive generations of the population, an impulsive control strategy known as pulsed chemotherapy for tumour control [4] and as pulsed application of pesticide or pulsed release of natural enemies for pest control [5]. Thus, the parameter θ is referred to as the dose timing response. The parameter q (i.e. the pesticide instant killing rate in Math.) is characterized by the pesticide efficacy of a particular dose applied at time point n + θ [4,5]. Thus, the parameter q, for convenience, is referred to as the dose response in the main text.
Note that the carrying capacity parameter K does not influence the dynamics of models EBHM and HRM [1]. Therefore, one of our main purposes is to reveal the effects of population dynamics (intrinsic growth rate r here), the dose response (parameter q here) and the dose timing response (parameter θ here) on the occurrence of paradoxical and hormetic effects.
Of course other factors such as environmental temperature, diet and pesticide or drug resistance, may also play important roles in inducing hormesis.
Through detailed theoretical analyses, we provide here the complete threeparameter spaces for the occurrence of the above phenomena, which provide guidance on experimental design and determine potential risks.

Dynamics of HRM with θ ∈ (0, 1)
The non-trivial equilibrium N * of model (S1.13) satisfies the following i.e., we have and denote the function F (N * ) as follows

The existence of equilibrium N *
In order to discuss the existence of the positive equilibria of model (S1.13), we only need to address the existence of positive roots of equation To do this, we denote with A = θ (θ−1)pe rθ−1 , please see more details from the last section of this SI for the definition of the Lambert W function. Obviously, we have A < 0 due to θ ∈ (0, 1) and the relations between A and −e −1 will play a key role in determining the sign of F ′ (N ) and the number of positive roots of the equation F (N ) = 0. Note that, according to the definition of the Lambert W function, A ≥ −e −1 ensures that N 1 and N 2 are well defined, and N 1 = Thus, at first we examine the equation Solving the above equation with respect to r, yields and consequently we denote the parameter sets Ω A>0 (r,p,θ) (or Ω A≥0 (r,p,θ) ) and Ω A<0 (r,p,θ) as the parameter spaces for A > −e −1 (or A ≥ −e −1 ) and A < −e −1 , respectively. Unless otherwise stated, similar notation has been used throughout.
Now we first discuss the signs of both functions F (N 1 ) and F (N 2 ), so for k = 0 and −1 we let Therefore, in order to ensure that W (k, A) are well defined for k = 0, −1, This means that r + ln(p) ≥ 4, i.e., we have pe r ≥ e 4 . For k = 0, we must have which holds true naturally for a ≤ −4. However, for k = −1, we must have which also holds true naturally for a ≤ −4.
Furthermore according to the definition of the Lambert W function we which depict the relations among three key parameters r, p and θ such that F (N 1 ) = 0 and F (N 2 ) = 0, respectively, denoted by Ω (r,p,θ) with r + ln(p) ≥ 4. Therefore, based on the equations given by (S2.5), the three-parameter space which determines the signs of the functions F (N 1 ) and F (N 2 ) can be implicitly given, as shown in Fig.S.1. This shows that the parameter spaces related to the signs of both functions are quite complex.
Furthermore, we have the following results.
] . (S2.6) According to the expression of N i (i = 1, 2), we have Therefore, g(N i ) < r(i = 1, 2). All these results confirm that This completes the proof of Lemma 1.

Proof.
If follows from the definition of the function F (N ) and the above notation that It follows from This completes the proof of Lemma 2.
Based on the above results, we first discuss the monotonicity of function  and −e −1 ≤ A < 0, then F ′ (N ) < 0 for all N ∈ Proof. Taking the derivative of F (N ) with respect to N , yields ))] (S2.10) To determine the sign of F ′ (N ), the signs of two parts, i.e., 1 − rθ h(N ) > 0, and consequently F ′ (N ) < 0. Therefore, we will focus on the case 1 − rθ K N < 0 (i.e., N > K rθ ) in the following. Letting F ′ (N ) = 0, we have h(N ) = l(N ), that is Note that l(N ) > 0 for all N > K rθ , and lim Moreover, by a simple calculation, we have , and it follows from It follows from the properties of the function xe Case b. −e −1 ≤ A < 0, i.e., for equation (S2.11) there exist two positive roots N 1 and N 2 .
It follows from This completes the proof.

Positive equilibrium and its stability
According to Lemma 2, the necessary condition for the existence of the positive equilibrium is pe r > 1 (i.e., F (0) > 0). Thus, we assume that pe r > 1 (i.e., r + ln(p) > 0 or a < 0) in the following. Based on the monotonicity of the function F (N ) the following cases for the existence of positive roots of the function F (N ) = 0 are discussed below.  Furthermore, if the positive roots exist, then we denote them as N * . The relations among N * i (i = 1, 2, 3) and K are quite important to show the effects of control measures on outbreaks of the pest population, and to depict the occurrences of a paradox and hormesis due to chemical applications.
According to the value of rθ, we consider the following several cases.
For this case we have the following relations: Based on the above inequalities, we can obtain that if F (N 1 ) > 0, then Note that if p = 0 and θ = 0, then F (K) = F (0). Thus, for this very special case the carrying capacity K is a unique positive equilibrium of the model. However, the values of positive equilibria, N * 1 , N * 2 , N * 3 , could be largely changed once a single pesticide application has been applied within the generation. The formula clearly shows how the dose timing (i.e., θ) and the dose response (i.e., p) of the pesticide applications affect the sign of F (K). To address this in more detail, we have the following results.
Lemma 3 Assume that r > 1 and p ∈ (p 2 c , p 1 c ), then there exists a unique threshold parameter θ c ∈ (0, 1) such that

Proof.
Letting and solving it with respect to θ, yields It is easy to see that θ c < 1. To ensure θ c ∈ (0, 1), we first consider the equation r(1−p)+ln(p) = 0 with respect to p. Rearranging r(1−p)+ln(p) = 0, yields the following equation Solving the above equation with respect to p and according to p ∈ [0, 1], yields two roots Similarly, according to the properties of the function −re −r and the Lambert W function, we see that p 2 c = 1 for r ∈ (0, 1] or r = −W (0, −re −r ) for r ∈ (0, 1]. Moreover, by a simple calculation for r ̸ = 1, we have dp 2 which indicates that dp 2 c dr < 0 for r > 1, and consequently, we have p 2 c < 1 for r > 1. Moreover, the function r(1 − p) + ln(p) with respect to p reaches its maximal value at p = 1 r with r > 1. These results confirm that the p 2 c is well defined for all given r > 1, and consequently, the threshold θ c given in the above is also well defined for all p ∈ (p 2 c , p 1 c ) and r > 1. This completes the proof.
Lemma 3 reveals that for a certain range of instant killing rate and a relatively large intrinsic growth rate r (i.e., r > 1 here), the different timings of pesticide applications may either successfully suppress the density of the pest population below the carrying capacity or cause pest population outbreaks and even resurgence with a density larger than the carrying capacity K, with the occurrence of hormetic and paradoxical effects.
Note that we can also solve r from the equation .
Thus, if we consider F (K) as a function of r, then we have F (K) > 0 for all r > r 4 c and F (K) < 0 for all r < r 4 c .
For this case we have the relations For this case, the relation between rθ and A 1 (or rθ and A 2 ) plays an important role. Note that A is a monotonically increasing function with respect to r, which indicates that A 1 is a monotonically decreasing function with respect to r, while A 2 is a monotonically increasing function with respect to r. Moreover, A 1 = A 2 = 2 at r = r c and + 2, which shows that r c θ = 2 when p = θ 1−θ . Furthermore, we have the following main results.
Proof If r c θ ≤ rθ ≤ A 1 ≤ 2, then we first consider the equation and solving the above equation with respect to r, yields ] .
It is easy to show that R c ≤ 0 (i.e., r c ≤ r ′ c ) and equality holds true if and only if θ (1−θ)p = 1 (i.e., p = θ 1−θ ). Similarly, we can show that the other results are true, and the proof is completed.
Note that there exists a unique r 4 c = − ln(p) (1−p)(1−θ) such that F (K) = 0 at r = r 4 c . Moreover, F (K) > 0 for all r > r 4 c and F (K) < 0 for all r < r 4 c . Further, we consider the relation between F (K) and F (K/(rθ)). where r 5 c and r 6 c are given in the following.
Proof It is easy to see that We define the function that is, It follows from F (K) > −1 that H(r) = 0.
Therefore, based on the above discussions, we conclude that the relationships among those threshold values r 2 c , r 3 c , r 4 c , r 5 c and r 6 c are quite complex, and Fig.S.4 shows one of the possible cases.
Based on the above inequalities, we easily have and the results can be discussed as those in case B1. That is, we can obtain In this case, we have the following relations: Thus, if F (N 1 ) > 0, then equation F (N ) = 0 has only one root N * 3 with N * 3 > N 2 > K; if F (0) > 0 and F (N 2 ) < 0, then equation F (N ) = 0 has only one root N * 1 with N * 1 < N 1 < K; if F (N 1 ) < 0 and F (K) > 0, then equation F (N ) = 0 has three roots N * 1 , N * 2 , N * 3 with N * 1 < N 1 < N * 2 < K < N 2 < N * 3 ; if F (N 2 ) > 0 and F (K) < 0, then equation F (N ) = 0 has three roots
By using methods similar to those used above, we consider the following two sub-cases.
Because rθ ≤ A 2 , therefore, N 2 ≥ K. Thus, we have Obviously, this case is the same as case B2.2 and can be similarly discussed.

The local stability of equilibrium N *
Based on the above subsection, there are at most three positive roots for equation (S2.9), denoted by N * 1 , N * 2 , N * 3 , respectively. Cases Note that for the first case rθ < rcθ we only need F (K) < 0 or F (K) > 0. rcθ 1 = max{1, rcθ}, Taking the derivative of f (N ) with respect to N , yields Thus, if pe r > 1, then N * 0 is unstable. This indicates that the pest population could be eradicated by spraying pesticides provided that the instant killing rate q is large enough (i.e., the survival rate p is small enough).
For the stability of positive equilibria N * , we have It follows from the monotonicity of the function F (N ) that we have F ′ (N ) > 0 for all N ∈ (N 1 , N 2 ), which indicates that F ′ (N * 2 ) > 0 due to N * 2 ∈ (N 1 , N 2 ). It follows from (S2.10) that g ′ (N * 2 ) > 0, and consequently we have f ′ (N * 2 ) > 1. All these results confirm that N * 2 is unstable if it exists. For the stabilities of N * 1 and N * 3 , we need |f ′ (N * )| < 1, i.e., −2 < N * g ′ (N * ) < 0. Note that if N * 1 (or N * 3 ) or both are the positive equilibria of model (S1.13), then it follows from the proof of Theorem 1 that the signs of both F ′ (N ) and g ′ (N ) are equivalent, and then we conclude that g ′ (N * 1 ) < 0 and g ′ (N * 3 ) < 0. Therefore, in order to show the stabilities of both N * 1 and N * 3 , we only need to show whether the inequality −2 < N * g ′ (N * ) holds true or not. Thus, we first consider the following function D(N ), i.e., that we have

Solving this equation for h(N ), we have
Substituting h(N ) into (S3.12), we can get Further, solving the equation D(N ) = 0 with respect to N , yields two roots as followsN Thus, based on the above discussion, we have D(N ) < 0 for all N ∈ (N 1 ,N 2 ), It follows from r + ln(p) > 0 (i.e., pe r > 1) thatN 2 > 0, and it is easy to see that r + ln(p) − 2 > 0 is equivalent toN 1 > 0. Therefore, based on this equivalence, for the local stability of all possible positive equilibria we consider the following two cases.
Case S 1 0 < r + ln p ≤ 2 In this case,N 1 ≤ 0 and D(N ) < 0 for all N ∈ (0,N 2 ). By a simple calculation, we have Let 2Z(r, p)), 2Z(r, p)), Proof We first consider the following equation N 1 =N 2 , i.e., Rearranging the above equation gives which indicates that Moreover, for all 0 < r + ln(p) ≤ 2, we have and this shows that equation (S3.13) is well defined according to the domain of the Lambert W function. Thus, we can obtain the following relation: from which the upper branch of W (0, A) of (S3.14) with respect to Z(r, p) is well defined and reveals the relationship among three parameters r, p, θ such that N 1 =N 2 . Therefore, r,p,θ and N 1 ≥N 2 ⇔ (r, p, θ) ∈ C Ω 1 r,p,θ . This completes the proof.

Theorem 2 For Case
In order to determine the stability of N * 3 , we need the relations between N 2 andN 2 . Let 2Z(r, p)), By employing the same methods as those in Lemma 7, we have the following results.
We list all possible cases for the stability of positive equilibria and their parameter spaces in Tab.S.2.
By a simple calculation, we have ] .
Now we discuss the sign of ∆ g θ . Denote .
According to the above discussion, we can conclude that (ii θ ) if Z < Z 1 or Z > Z 2 , there are two rootsN 5 andN 6 of equation then g θ (N * ) > 0, wherē , .
Because Z 2 > 1, therefore,N 6 >N 4 , when Z ≥ Z 2 . Due to Z 1 < 1, thus, This completes the proof. In summary, we have the following main results.
The effects of parameter θ on the positive equilibria, hormetic and paradoxical effects are listed in Tab.S.3.

Multiple control actions within each generation
Without loss of generality, we assume that the chemical control tactics have been applied m times within each generation [n, n + 1], i.e., there exist θ i (i = 1, 2, · · · , m) with n ≤ n + θ 1 ≤ n + θ 2 ≤ · · · ≤ n + θ m ≤ n + 1, such that the control measures have been applied at n + θ i with a proportion q i of the pest being killed. Thus, the survival rate after each chemical control is p i = 1 − q i . For convenience, we denote θ 0 = 0, θ m+1 = 1 and p 0 = 1.
Therefore, by employing the same methods proposed before we can extend the Beverton-Holt model with multiple controls within each generation as In particular, if m = 1 then we have The Ricker model with multiple controls within each generation can be generalized as follows and In particular, we have ) and A 4 = exp ( rθ 4 − r K N n (θ 1 + (θ 2 − θ 1 )p 1 A 1 + (θ 3 − θ 2 )p 1 p 2 A 2 + (θ 4 − θ 3 )p 1 p 2 p 3 A 3 ) ) .

Some important definitions related to the Lambert W function
Definition The Lambert W function is defined to be a multivalued inverse of the function z → ze z satisfying LambertW (z) exp(LambertW (z)) = z.  [4] Panetta, J.C. A mathematical model of periodically pulsed chemothera-