BACTERIA–PHAGOCYTE DYNAMICS, AXIOMATIC MODELLING AND MASS-ACTION KINETICS

Axiomatic modeling is ensued to provide a family of models that describe bacterial growth in the presence of phagocytes, or, more generally, prey dynamics in a large spatially homogenous eco-system. A classification of the possible bifurcation diagrams that arise in such models is presented. It is shown that other commonly used models that do not belong to this class may miss important features that are associated with the limited growth curve of the bacteria (prey) and the saturation associated with the phagocytosis (predator kill) term. Notably, these features appear at relatively low concentrations, much below the saturation range. Finally, combining this model with a model of neutrophil dynamics in the blood after chemotherapy treatments we obtain new insights regarding the development of infections under neutropenic conditions.

1. Introduction.The continuous fight between bacteria and the white blood cells that belong to the innate immune system occurs in a large class of animal species (from insects to humans 1 [1]).Indeed, this process is of tremendous clinical importance as these white blood cells provide the main defence against bacterial infections in the acute stage [2,3].Yet, there are very few mathematical models that describe its basic properties.As a first step for exploring the very involved in-vivo fight, where normally the phagocyte 2 concentration is strongly coupled to the bacterial load, one needs to understand the in-vitro dynamics where the phagocyte concentration is fixed.In [6] a model for the bacterial growth in the presence of a constant concentration of neutrophils (white blood cells that eliminate bacteria by phagocytosis) was proposed, and its clinical implications were explored.We focus on neutrophils as they are the most abundant type of phagocytes (about 70% of the white cell blood count in an adult) and hereafter phagocytes and neutrophils are used interchangeably.Notably, beyond being a basic building block to other in-vivo models, this model may be relevant to some specific medically important conditions under which it is reasonable to assume that the phagocyte concentration remains constant.These conditions are characterised by either insufficient supply of phagocytic cells by the bone-marrow or by ample supply of malfunctioning phagocytic cells (for details, see [3,6] and the references therein).Such conditions are associated with high risk of infection [7,8].Another natural setting for considering constant phagocyte models are in-vitro tests that are performed to detect phagocyte malfunction (e.g., [9]).Therefore, understanding the bacteria-phagocyte interaction at constant phagocyte concentration is important clinically.Biologically, similar conclusions apply to population modeling of prey dynamics in a well mixed environment with approximately constant predator concentration.Mathematically, we view this very simple model as a carefully constructed building block to other models.Indeed, the bacteria-phagocyte interaction is a basic ingredient in the complex innate immune response to bacterial infection (see e.g.[10,11,12]).
Medically, the blood neutrophil concentration is viewed as a defining measure for assessing the risk of infection.Patients having concentrations below the threshold value of 5 × 10 5 neutrophils/mL in the blood (called neutropenic patients) exhibit high risk of infection, with a sharply increased risk associated with both the duration and the severity of the neutropenia [7,8].Mathematically one may argue with the use of the term "threshold" here, as there is no single critical value below which all patients become sick and above which all remain healthy (and we will later on see that there is a mathematical reasoning that explains this observed variability).Yet, we use this term here as it reflects the common view by the medical community: severe neutropenia is medically defined by this single threshold value.Recently, it was proposed in [13,14] that in-vitro bacteria-phagocyte experiments explain this threshold effect.The authors conclude from their experiments that a single critical neutrophils concentration is required to predict the final outcome of the bacteria-neutrophil interaction, regardless of the initial bacteria concentration.Moreover, they showed that their fitted in-vitro critical neutrophil concentration is comparable to the medically defined one.Here, to better explain the meaning of this in-vivo threshold level we construct a model that presents a simplified view of the in-vivo dynamics under neutropenic conditions.To this aim, we introduce into the model [6] a typical externally controlled variation of the neutrophils dynamics after chemotherapeutical treatments as described in [15].By combining these two models we gain some new insights into the phenomenon of sharp transition to fast development of acute infection as the neutrophil levels are decreased.
In [6] we argued that by modeling the bacteria-neutrophil interaction using the axiomatic approach a more adequate model to the bacteria-phagocyte interaction at constant phagocyte concentration may be constructed.Here, we provide the mathematical formulation and justification of the main claim made in [6].Under reasonable biological assumptions on such models there are exactly two simplest possible robust bifurcation diagrams (monostable and bistable dynamics).The model construction in the most general form that is consistent with the biological theory, so that its predictions are not form-and parameter-specific.Notably, predictions are viewed here in the broad sense of qualitative analysis of dynamical systems.Concretely, as pointed out in general biological terms in [6] and mathematically defined and proved here, we are able to classify all the possible robust sequences of bifurcations that occur when the key parameter, the phagocytic cell concentration, is varied.Surprisingly, we show that at low bacteria concentrations there are exactly two such sequences: forward and backward transcritical bifurcations.This result provides new predictions that are somewhat counter-intuitive for experimental biologists working in the field.
This axiomatic model construction fits within the general dynamical-system mathematical framework of obtaining robust qualitative predictions regarding biological phenomena by constructing families of models with certain properties [16,17,18,19,20].For example, symmetries of the interaction terms (group theoretic constraints) have been often used to predict dynamical properties of biological systems, see e.g.[21].Similarly, the class of monotone systems [22] and its generalization e.g, [23,24,25] were utilized to analyze some models of biological processes and of reaction networks [26].Other qualitative properties, from the existence of robust homoclinic cycles [27] to characterizations of partial derivatives [28] have been utilized in a similar manner.This modeling approach has been successfully applied to specific medically-relevant biological problems that are related to the innate immune system [15,29].
The axiomatic model construction begins in Section 2 with the selection of the relevant dependent variables (here the bacteria concentration, all other key ingredients including phagocytes, nutrients, opsonins, cytokines etc are assumed to be of constant concentration), independent variables (here time, spatial effects are neglected) and the evolution law type (here, deterministic ordinary differential equation, neglecting, for example the bacteria and phagocyte age distribution, stochastic effects, memory effects, delays etc.).We postulate that these assumptions are adequate for the relevant in-vitro setting, as in [13,14].Then the model is formulated in terms of the type of the biological interactions that govern the evolution of the dependent variables (here, the birth, influx, and elimination of bacteria) and the dependence of these interactions on the variables and on the key parameters.To ensure generality, these dependencies are formulated in terms of monotonicity properties of the various interaction terms and by assuming a form of these terms for sufficiently large and small values of the variables.
In Section 3 we analyse the dynamical properties of the resulting class of models.We classify all possible types of bifurcation diagrams and establish that at small bacterial concentrations only two types of behaviour are possible.In Section 4 several explicit models are presented.Some important properties of the model of [6] are listed and this model is parameterized by using the data from [14].We then discuss two other classes of commonly used models for prey-predator dynamics that do not satisfy some of our axioms.Section 5 describes a new model, in which we combine the bacteria-neutrophil in-vitro model with a previously published model of the neutrophils blood dynamics after chemotherapy [15].Finally, Section 6 presents the summary of our findings and discussion of their implications.

2.
The model axiomatic construction.We postulate that the rate of change of the bacteria concentration b = [#bacteria/ml], (b ∈ Ω ⊂ R + ) in a well mixed chemostat may be described for a sufficiently large population by an ODE model: where the terms B, D represent the "birth" and "death" processes3 that depend on two parameters, the predator/killer concentration n and the bacterial flux4 s into the region.Here n may represent the neutrophils concentration as in [6], or concentration of other entity that destroys the bacteria (such as antibiotics and antimicrobial peptides [30]).The first two terms f (b; s) := B(b; s) − D(b) correspond to the net rate of change of bacteria in the absence of neutrophils (natural bacteria dynamics) under the given chemostat conditions.The third term D K (b; n) corresponds to the kill term of the bacteria in a unit volume per unit time.Chemostat means here an ideal dish in which the bacteria environment remains unchanged forever.In reality, in-vitro bacteria-phagocyte interactions are tested in small vibrating tubes at body temperature and last about 60-90 minutes -a time scale at which the bacteria environment is almost unchanged and the bacteria can maintain exponential growth [31].At longer time scales it is known that the phagocytic function is reduced and the assumption of constant environmental conditions is not applicable to the common experimental setting.Nonetheless, we propose that the above model is adequate for describing these experiments that are stopped before a significant change in the environment is observed.

Empirical Properties of the Rate Functions.
There are several robust biological observations regarding the growth of bacteria: the natural (no phagocytes) bacteria growth model always leads to a limited growth curve.Away from the maximum capacity regime the growth curve is monotone and is well-described by one dimensional dynamics [31].The kill term typically monotonically increases with the bacteria and the phagocyte concentrations up to some saturated value.These robust observations imply specific mathematical properties for the birth and death terms of the bacteria model that we summarize by Assumption A1-A6: Indeed, since the time scale of the current model is comparable to the time scale associated with the division and the natural death/kill processes, the birth and death rates per bacteria must be finite.
A 2. When there is no external source of bacteria and no bacteria are present in the body, there is no spontaneous growth or natural death of bacteria Similarly, killing occurs only when both the bacteria and the neutrophils are present A 3. The natural bacteria dynamics obey a limited growth model [31]; The birth term is monotonically increasing with b, s The net bacterial growth rate changes sign at the maximal capacity concentration b f (s) Standard growth models such as the logistic or Gompertz models [31] satisfy A3.
A 4. The bacteria killing term is monotonically increasing with b and n [32,33]: and the local killing rate (i.e. the killing term slope) is monotonically increasing with n: This monotonicity assumption is expected to be widely valid in the predator-prey context up to moderate concentrations, at which, for example, spatial effects may become significant and may lead to non-monotonic effects on the population level.
A 5. A sufficiently large neutrophil level can overcome limited bacterial challenge [32,33,34,35].Namely, there exist positive b max , s max values such that for all (b l , s l ) ∈ (0, b max ) × (0, s max ) there exists n 0 (b l , s l ) > 0 such that for all n > n 0 (b l , s l ) the bacteria die whenever its concentration and influx are below b l and s l respectively This assumption reflects the clinically and experimentally established efficiency of the phagocytes in overcoming limited bacterial challenge.In section 4 we provide concrete examples for interaction terms that satisfy these assumptions.A 6. Apart from the specific properties listed in A2-A5, the functions (B(b, s) − D(b)), D K (b, n) are in "general position".In particular, the following conditions are satisfied This last assumption states that there are no additional constraints apart from A1-A5 and hence that there should be no accidental cancellations between higher order derivatives.Put differently, assumptions A1-A5 mean that we may write where the hatted functions satisfy various inequality conditions.Assumption A6 says that for finite arguments the hatted functions are "generic" smooth functions, namely, they do not satisfy any additional specific constraints.Thus, we conclude that the class of the corresponding hatted functions is a residual set, i.e. it is the complement of a meager set in the sense of Baire category [36].
Even though the non-degeneracy assumptions that are listed in A6 arise from mathematical reasoning, we do expect that concrete models that adequately reflect the biology will satisfy these conditions.It simply reflects the belief that degeneracies should not be created incidently, and in any case can be removed by small changes in the model.
3. Dynamical Properties.The main mathematical results of the manuscript are formulated below.Roughly, these imply that models of the form of equation ( 1) have exactly two possible behaviors at small b: a single transcritical bifurcation either supercritical or subcritical (sometimes called forward or backward bifurcation [37,38,39]) and quite simple behavior at higher b values (see Appendix for additional insights).
Theorem 3.1.The zero flux bacterial growth model, namely equation (1) with F (b; n, s = 0) satisfying assumptions A1-A6, has a transcritical bifurcation at b = 0, at some unique value n = n trc > 0. All other bifurcation points (if exist) are saddlenode bifurcation points, and occur at b > 0.Moreover, 1. ∂ 2 F ∂b 2 (0; n trc , 0) < 0 implies that the transcritical bifurcation is supercritical (forward) in (n trc − n) and that there are additional even number (possibly zero) of saddle-node bifurcation points.2. ∂ 2 F ∂b 2 (0; n trc , 0) > 0 implies that the transcritical bifurcation is subcritical (backward) in (n trc − n) and that there are additional odd number of saddle-node bifurcation points (in particular at least one such bifurcation point exists).
See the Appendix for details of the proof.The main idea of the proof is that under assumptions A1-A6, the nontrivial equilibria set may be represented as a graph n = n * (b) that intersects the b axis at b = b f and the n-axis at n = n tcr .The sign of ∂ 2 F ∂b 2 (0; n trc , 0) determines the slope of this curve at the intersection with the n axis.The inequality (12) implies that this slope also depends on the kill term parameters.Finally, the equilibrium point (b, n * (b)) is stable (respectively unstable) when n * (b) is decreasing (respectively increasing) with b, and the segment emanating from n = 0 at b = b f is always stable and decreasing.
When the bacteria influx s is positive and small, by assumptions A3-A4 and by utilizing the Implicit Function Theorem (IFT), the stable and unstable equilibrium branches of F (b; n, 0) may be continued in s to the corresponding branches of F (b; n, s) with the same stability type.Moreover, the direction of the shift in the position of these branches with s is determined by their stability: Lemma 3.2.A zero flux stable (respectively unstable) equilibrium point (b * (0); n * ) is pushed upward: for a fixed n * , for sufficiently small influx s, b * (s) > b * (0) (respectively downward, so b * (s) < b * (0)).Similarly, the equilibrium point curve is pushed to the right: for a fixed b * belonging to the equilibria curve, for sufficiently small influx s, n * (s) > n * (0).
Proof.: Notice that by implicit differentiation of F (b * (s); n * , s) of (1) at a fixed n * satisfying F (b * (0); n * , 0) = 0, we obtain dF (b * (s); n * , s) = ∂ s F ds + ∂ b F db. Thus, when (b * (0); n * , 0) is not a bifurcation point of the zero-flux model (so ∂ b F (b * (0); n * , 0) = 0 ) we get by the implicit function theorem that for small s, ds is positive at a stable fixed point (where ∂ b F (b, n, s) < 0) and similarly it is negative at an unstable fixed point.Thus for a non-bifurcation point, for small s, b * (s) is a See Figure 2 for the resulting two types of bifurcation diagrams that appear at small influx s.Notice that the trivial zero-flux solution b = 0 is pushed by the bacterial influx to negative values for n < n trc and to positive values for n > n trc .Thus, in the subcritical case, the transcritical bifurcation splits to two saddle nodes, only one of which appears at the biologically relevant positive quadrant regime.Theorem 3.3.For sufficiently small bacterial influx s > 0, the bacterial growth model (1) satisfying assumptions A1-A6 has an even number of saddle-node bifurcation points and no other bifurcations in the positive quadrant.If ∂ 2 F ∂b 2 (0; n trc , 0) > 0 then there are at least two saddle-node bifurcation points.
See Appendix for details of the proof.This theorem essentially follows from Theorem 1, the Implicit Function Theorem, Lemma 3.2, and the observation regarding The above two theorems lead us to the following main conclusion.There are exactly two "simplest" type of behaviors, where "simplest" models are those that have the least number of bifurcation points: Corollary 1.The simplest model of the form (1) satisfying A1-A6, with s = 0 (respectively s > 0) has exactly two possible bifurcation diagrams: 1b (respectively Fig. 2b).The dependence on s of the bifurcation point near (b, n) = (0, n trc ) has the typical square-root singularity that arises when the transcritical bifurcation symmetry is broken.Indeed, the saddle-node bifurcation point of the perturbed transcritical normal form 6 . Thus, small changes in the bacteria influx s change the location of the bifurcation point dramatically (since 4. Several explicit models.We present three models -the first is a model that was introduced in [6].It is a valid model (hereafter, a model that satisfies assumptions A1-A6) that is derived phenomenologically and presents the above simplest form of dynamics.The other two models correspond to commonly used models that are not valid (violate A5 and A6).Interestingly, these models indeed miss biologically important behaviors.The first non-valid model utilizes a quadratic "mass-action" kill term [13,14].We show that when the bacteria natural growth rate function is monotonically decreasing in b such a model cannot exhibit bistability, a phenomenon that is experimentally observable and has tremendous clinical and biological implications [6].The other non-valid class of models are the ratiodependent models.Such models miss the important biological effect of the existence of a minimal phagocyte concentration threshold [13,14,40].

4.1.
A valid phenomenological model.In [6] the following model was introduced and analyzed where the first two terms are the natural bacterial growth while the last term is the killing term, see e.g., [41,42,17] for a discussion of models with similar algebraic structure).It is easy to show that this model satisfies assumptions A1-A6 for sufficiently small bacterial influx ( 0 ≤ s ρ/β) provided the six parameters (ρ, δ, β, α, γ, η) are all positive and satisfy ρ > δ (to get positive maximal capacity at The first three parameters govern the natural dynamics of the bacteria: ρ and δ control the natural linear growth/death rates of the bacteria and β controls the natural saturation of the bacterial growth rate at high concentrations.The other three parameters control the killing of the bacteria by the neutrophils: α is the neutrophils' bacterial killing rate at low concentrations, and γ and η control the saturation and interference in the killing rate as the concentrations of b and n are increased (see A4).All of these parameters depend on both the bacterial strain and the environmental setting.For example, in an experimental set-up, the serum content and the form of the toxic clearance affect them.In fact, these parameters represent global effective responses to a large complex network of molecular and cellular processes that are associated with the bacteria natural growth and with the capture and killing of the bacteria by neutrophils.A calibration of these parameters with extensive experimental data and with the corresponding molecular markers may be important for identifying the main effect of these cellular processes on the bacteria-neutrophil population dynamics in health and in disease.
This model exhibits the simplest type of dynamics as listed in Corollary 1.At s = 0, for the parameter regime it is monostable whereas when the opposite inequality occur it is bistable [6] (to establish the above condition, one finds when ∂ 2 F (b,n,s) | (0,ntrc,0) < 0, namely when the transcritical bifurcation point is non-linearly stable).Lemma 3.2 implies that for sufficiently small non-zero flux s the same conclusion holds.
The inequality demonstrates that if one fixes the parameters (ρ, δ, η, α) to the biological regime (so the right hand side of the inequality is positive), then, for sufficiently small 7 β and/or for sufficiently large γ bi-stability emerges.Smaller β corresponds to saturation at higher levels in the bacteria growth, namely, a situation by which there is more space and nutrients supply to the bacteria.Larger γ corresponds to saturation at lower bacteria levels in the kill term, meaning that the neutrophils reach their maximal killing capacity at a lower bacterial concentrations.The inequality shows that more favorable growth conditions for the bacteria and/or less efficient phagocytosis at high bacterial concentrations lead to bi-stability.Indeed, Figure 1 and 2 compare the bifurcation diagrams of (13) at two parameter sets that differ only in their β values: when β = 9•10 −9 (Figures 1a,2a) we see monostability whereas at the lower β = 10 −9 (Figures 1b,2b) bi-stability emerges.In fact, Figure 3 of [6] shows that for relevant parameter regimes the bistable scenario is expected to be more common than the mono-stable one.
Indeed, we propose that the experimental results of the in vitro S. epidermidisneutrophil interaction in [14] are consistent with the bistable behavior.Since the growth rate data for this experiment are not available to us, we fix the natural linear bacterial growth rate to be as in [14] and choose the elimination rate to be considerably slower (by two orders of magnitude) and the growth rate saturation to be 10 −8 [1/#bac/ml] (namely, we assume that the bacterial population capacity is 10 10 [#bac/ml] see e.g., [31,43]).We then fit solutions of the ordinary differential equation to the data by adjusting the killing term parameters of (13) using a least square cost function with multiplicative measurement model [44] and the Matlab implementation of the multidimensional downhill simplex method for the minimization (see e.g., [45,46]).We find that the fitted kill-term parameters are well within the bistable regime (and this conclusion is independent of the particular choice for the natural bacterial growth parameters as long as β ≤ 10 −8 ).For example, for these parameters, condition (14) implies that to achieve monostability γ needs to be decreased by a factor larger than 8.In fact, the 95% confidence interval for the estimated γ, evaluated using the bootstrap percentile method [47] with 1000 samples, lies entirely in the bistable region.This is true for the entire cube defined by the confidence intervals of α, η, and γ.These conclusions hold for β ≤ 10 −8 .They fail for higher values of β, namely for population capacity smaller than 10 10 [#bac/ml].Figure 3 shows the bifurcation diagram for the fitted parameters at s = 0.01, 100, 400, 1500 [#bac/ml min] (see caption for the parameter values).
Let us note several important properties of this model at the positive small influx case that will be used in section 5. First we recall that by Theorem 3.3, for the subcritical (backward) bifurcaion case, the bifurcation point (b, n) = (0, n trc ) becomes a saddle-node bifurcation point.Asymptotic expansion for this point yields: Notice that y is positive for bistable parameter sets (see Eq. ( 14)).In particular, for the fitted parameter values, y 1.13 • 10 9 , so n sn (s) ≈ n trc + 5 • 10 3 √ s and the critical bacterial load at this point rapidly increases as √ 10 9 s. Figure 3 shows that for the fitted model the asymptotic expansion accurately predicts the bifurcation point location up to s = 400 [#bac/ml/min], whereas for s > 500 the higher order terms contribute.At s = 1500 [#bac/ml/min] there is a marked difference between n sn (s) and n trc .Notice that for a given n value there is a critical bacterial influx s c (n), such that for s > s c the neutrophil level n is beneath the threshold (n < n sn (s)) and thus infections develop to the acute stage.Using the above asymptotic expansion we find the asymptotic expansion for s c (n) (valid for n n trc ) For the parameters corresponding to [14], and for n n trc = O(1), the error terms are of order O(10 −9 s).These findings are important for understanding the development of infections under neutropenic conditions, see section 5.

4.2.
Models with simple mass-action kill term are not valid.Many predatorprey models use as a building block a quadratic kill term of the form D K (b, n) = αbn ("mass-action kinetics" assumption) where the natural bacteria dynamics is taken to be either linear or of limited growth.Such a kill term, which is linear in b, is considered to be a reasonable approximation to mimic competition and/or predation between species when the concentrations are much smaller than the saturation regime.For example, the classical SIR models [48] employ such quadratic interaction terms.While it is well established by now that saturation effects are important for large concentrations of populations [32,13,14], we show next that even for small concentrations these terms play an important role in preserving the robustness of such models.
Indeed, notice that with a quadratic kill term, the nontrivial fixed point may be simply found from In particular, if the bacteria natural growth rate function R(b; s) is monotonically decreasing with b (as expected, as is the case for the logistic and other common models for the bacteria natural growth) then so does n * (b).By the same arguments as in Lemma 6.4 we conclude that this branch must be stable.We conclude that when the natural bacterial growth rate is monotonically decreasing with b, a mass action kill term always results in a single monotonically decreasing stable branch of equilibrium.In particular, bistability is not possible.Such models violate only one of our assumptions in the axiomatic formulation: the inequality (12) of the "general position" assumption A6 is not satisfied.Thus, models with quadratic kill term correspond to a restricted class of models satisfying assumptions A1-A5, and for this restricted class, when R(b; s) is monotonically decreasing with b, only monostability is possible.
In experiments it is often further assumed that the bacteria are in a phase of exponential-growth, which means that the natural bacteria growth term is well approximated by a linear term B(b, 0) − D(b) = rb [13,14].Combining this linear growth with a mass action kill term (so db/dt = (r − nα)b) results in a degenerate bifurcation diagram.Indeed, this bifurcation diagram consists of a single vertical line emanating from n c = r/α (where the right hand side identically vanishes), and a single horizontal line -the n axis-that corresponds to the trivial solution b = 0.Such a diagram does not appear in textbooks and rightly so.It corresponds to a linear truncation of a model with a bifurcation, where we know that higher order terms must be included to determine the stability near the bifurcation point.In particular, the addition of a non-linear term of the form αnb 2 changes the qualitative behavior of the bacteria dynamics near n trc for arbitrarily small b (see Lemma 6.5 of Appendix 1).Put differently, near a bifurcation point non-linear terms that contribute to the normal form must be included.Neglecting the saturation in the kill term amounts to neglecting such non-linear corrections, and these corrections are the leading order terms at the bifurcation point.
We have thus shown that approximating the kill term of, for example, (13), by its leading order behavior, the quadratic term αbn, may lead to erroneous results near the bifurcation point (0, n trc ).In particular, we showed that if the bacteria natural growth rate function R(b; s) is either constant or monotonically decreasing with b, models with mass-action kill term are not bistable.Finally, we pointed out that the over-simplified model that combines a linear bacterial growth and a mass action kill term is highly degenerate violating both A3 and A6.4.3.Ratio-dependent models are non-valid.In the experimental practice of phagocyte-bacteria interactions only ratios were considered to be important [32,34,33,35].A ratio-dependent model for the kill term assumes that the bacteria population depends only on the initial ratio between the bacteria and the neutrophils concentration, namely, it consists of models of the form (here we take s = 0) When the kill term ) k is ratio dependent whereas the kill term of ( 13) is not ratio dependent (see e.g., [16,49] for predator-prey model with such ratio-dependent kill terms).
We first note that the ratio-dependent models have a non smooth (not C 2 ) kill term at n = b = 0, and thus violate assumption A1.Indeed, at n = 0 the kill term vanishes, and the origin is unstable (since R(0) > 0).On the other hand, for all n > 0, provided R(0) = K(0), the origin is hyperbolic with an attraction/repulsion rate (R(0)−K(0)).Thus, if K(0) = 0 these models are not smooth (taking K(0) = 0 solves this smoothness problem but leaves the origin unstable for all n, violating A5).
We also notice the following peculiarity of such models.When8 R(0) < K(0) A5 is formally satisfied.Yet, we notice that the origin attraction rate is independent of n, a property which defies biological intuition by which the bacteria eradication is faster as n is increased.In fact, as Theorem 3.1 shows, combining assumptions A1 and A5 imply that any valid model must have a positive threshold neutrophil concentration below which the origin is unstable and above which it is stable.As a consequence of the non-smoothness of the ratio-dependent models, these models do not have such a positive threshold neutrophil concentration.
The above considerations show that ratio-dependent models are inadequate for describing the behavior at small bacterial concentrations.

5.
A model of the in-vivo dynamics in neutropenic patients.To better explain the potential medical implication of the bistable behavior, we consider a clinical situation by which the neutrophil dynamics is externally controlled.In particular, we examine when an acute infection develops as a result of a neutropenic "cycle" -an event of temporary reduction in the neutrophil level.To this aim, we use the parameterised model ( 13) (fitted to the in-vitro data from [14] of S. epidermidis-neutrophil interaction experiments as in Fig. 3) and combine it with a model that faithfully describes the neutrophils blood levels, N , of patients undergoing chemotherapy with supportive treatments of Granulocyte Colony Stimulating Factor (G-CSF) injections [29,15].Our model for the in-vivo dynamics is defined by equating the neutrophils level in the tissue to their level in the blood and by taking a positive constant9 bacterial influx s.
Before continuing in describing the model and its predictions, we should stress that this combined model does not mimic the true in-vivo dynamics even under neutropenic conditions.Instead, we argue below that usually, starting a couple of days after the onset of the neutropenia, the neutrophils concentration in the tissue is smaller than their blood concentration.Therefore, this model provides an over optimistic scenario with regards to the development of infection for such patients.While this model is not realistic, it does provide important information -it shows the most optimistic plausible outcome of the infection development for a given prolonged neutropenia profile and thus highlights the key parameters that determine when acute infection definitely develops.
The model for the G-CSF and the neutrophils blood levels was derived and analysed in [29,15] [15] values, all the other parameters are as in [15] (notice that we divide the daily rates of table 2 of the supplement of [15] by 1440 to obtain rates per minute).N * = 5 • 10 6 cells/ml is a fixed scaling parameter.
it is combined with equation ( 13) to form the in vivo model The neutrophils (N ), as discussed extensively in this paper, are the main defense of the body against bacterial and fungal infections.G-CSF (G) is naturally produced by the body and is also available for injections as a supporting medical treatment.
It is known to control the neutrophils production in the bone marrow and the neutrophils delivery into the blood.It is commonly administrated to neutropenic patients by either 10 daily standard subcutaneous G-CSF injections (Filgrastim 5 μg/kg/day), by a single pegilated -pegG injection (Pegfilgrastim 100 g/kg), or, in special cases by a continuous infusion (Filgrastim 10 μg/kg/day).In [15] it was suggested that bi-daily G-CSF injections, that are commonly used in blood harvesting protocols [50], may be also beneficial to a certain class of patients.The functions inG(t) (G-CSF injections) and BN F (t) (neutrophils influx from the bone marrow) in this model are external forcing functions that are induced by the medical treatments.In [15] it was demonstrated that the following functional forms for these forcing terms adequately represent several different clinical data sets that correspond to different chemotherapy treatment protocols with N gshots subcutaneous G-CSF injections, given at a T f period (here, either daily or bi-daily, i.e., N gshots is either 10 or 20) The term inG(t) corresponds to the pharmokinetics absorption of the G-CSF injections into the blood, so with no treatment and for t < T gstart , inG(t) = 0.The term BN F (t) mimics the influence of the chemotherapy and the G-CSF treatment protocol on the bone marrow capacity to produce and deliver neutrophils to the blood.Normally, with no treatments, BN F (t) = B N is a positive constant that reflects the normally constant influx of neutrophils from the bone-marrow into the blood.This term is multiplied in the GN model ( 19) by an increasing saturable function of G, reflecting the observation that under normal conditions this influx is enhanced when the G-CSF blood concentrations are high (as a result, for example, of an infection).The chemotherapy damages the bone marrow and thus makes the basic influx of the neutrophils into the blood plummet to low values (B nadir ) for a few days till its recovery.Fitting the form (20) of the BN F (t) to several clinical data sets shows that its parameters are treatment specific as they reflect the severity and character of the chemotherapy protocol [15] (one may suspect that these are also patient specific, but so far there are no clinical studies examining the inter-patient variability of these parameters).
Here, we employ a realistic set of parameters (see Table 1) of a patient undergoing chemotherapy treatment and experiencing as a result of it a prolonged severe neutropenia.The control parameter AM C = k N EF • B N adir N * •D N denotes the acute marrow capacity of the patient.This non-dimensional parameter controls the neutropenic response to G-CSF treatments, and with maximal G-CSF adminstration, N ≈ N * AM C. Therefore, AM C values above 0.2 correspond to cases by which G-CSF treatments suffice to keep the patient away from the neutropenic regime, while AM C < 0.1 corresponds to cases in which even with intensive G-CSF treatments the patient will suffer from sever neutropenia (see [15] for details).Here the AM C is set to 0.2.With no G-CSF support, this patient would suffer from approximately 10 days of severe neutropenia whereas standard daily injections (a standard injection consists of a dose of 5 μg/kg of Filgrastim) or bi-daily standard injections (a therapy advocated in [15] for patients at risk) both lift this patient's neutrophil counts to above the critical level.Figure 4a-c shows the bacteria, neutrophil and G-CSF dynamics of this patient with no G-CSF support (cyan curves), with 10 daily standard G-CSF injections (blue curves) and with bi-daily standard injections (green curves).A rough estimate for the achieved neutrophil levels at the nadir may be found by noticing that the injections bring the G − CSF value to around 30, 000pg/ml, so, by setting this value of G in the second equation of (19) and plugging in the parameter values of Table 1, the quasi-steady state neutrophils levels become N nadir ≈ 0.87 • N * • AM C ≈ 8.7 × 10 5 .Figure 4b shows that with the bidaily injections the neutrophil levels at the nadir indeed stabilizes to around N nadir (oscillates between (8.5−8.8)×10 5 neutrophils/ml) whereas with the daily injections the neutrophil levels are lower, oscillating between (5.7 − 7.7) × 10 5 neutrophils/ml.The red line in Figure 4b denotes the commonly used medical threshold level of N c = 500 × 10 3 #neutrophils/ml, demonstrating that according to current medical wisdom the standard G-CSF therapy of daily injections should suffice to abrogate the severe neutropenia for this patient.On the other hand, both the daily and bi-daily injections leave the patient below the 1, 000 × 10 3 #neutrophils/ml line for about 10 days, and according to [51], with such an extended period of neutropenia, patients have 20% risk to develop an infection.Such statistical statements imply that there exists inter-patient variability and that, as alluded in the introduction, the value N c is not really a strict universal threshold value.In (d) these trajectories are projected onto the neutrophils-bacteria plane.In (e) the response of the patient receiving bi-daily injections to an increased bacterial influx is demonstrated.When the influx is increased from s = 100 (green, as in (a)) to s = 400 (red) acute infection develops.When the kill term parameter α is increased by 20%, a bacterial influx as large as s = 3500 [#bac/ml/ min] results in an infection that is under control (magenta curve).
Here, we propose that the bistable model of the bacteria-neutrophil interaction may explain this variability and supply a more accurate, patient and infection specific interpretation of the threshold effect.Indeed, we assert that the value of n sn (s) ≈ n trc + O( √ s), see equation ( 15), provides a good working threshold value for the onset of infection.If, for all time, n(t) > n sn (s) and b(0) < b sn (s), then the combined model predicts that no acute infection will develop.On the other hand, when n(t) drops below n sn (s) the bacteria concentration increases, and, if it crosses b sn (s) the infection will usually increase to the acute phase (defined here as b > 10 8 , the horizontal black lines in Figure 4 (a) and (e) ).In particular, if n(t) drops below n trc for some non-negligible time10 even the tiniest bacterial influx will cause acute infection.
To demonstrate these claims, Figure 4a shows the resulting bacteria dynamics for the three G-CSF treatment profiles with a small bacterial influx of s = 100 [#bac/minute].We see that to overcome an infection of the S. epidermidis bacteria with the neutrophils having a kill-rate as in the experiments of [14], the patient needs to receive bi-daily injections, while the standard daily injections are insufficient.Indeed, we find that with these parameter values and with the influx parameter s = 100 [#bac/ml/min], the neutrophil threshold parameter is n sn (100) ≈ 8.4 × 10 5 [#neutrophils/ml], and since both the untreated patient and the patient receiving daily injections have substantial periods at which their neutrophil counts are below this value and in fact even below n trc , both will definitely develop acute infections.We see that here the bi-daily injections keep the patient very close to the border of instability, meaning that to provide a better assessment of the infection outcome for this patient a more realistic model is needed.Figure 4d shows the phase-space plots of these three profiles with the corresponding bifurcation curve (black dashed line).
Increasing the bacterial influx for this patient to s = 400 [#bac/ml/min] increases the neutrophil threshold value to n sn (400) ≈ 8.8 × 10 5 [#neutrophils/ml, and the bi-daily injections are indeed no longer sufficient to overcome such an invasion of bacteria -the red curve in Fig 4c shows that the bacteria population increases to the acute stage with this larger bacterial influx.
We conclude that for a given patient profile (the profile consists of the neutrophil dynamics with treatment and the neutrophils kill rate parameters) there exists some critical value of bacterial influx s c (n(t)) above which acute bacterial infection definitely ensues.This value of s c provides an important information: potentially (after more clinical data is gathered regarding the common values of s), it can tell us whether the patient is able to overcome the internal bacterial population in the body and whether it can withstand contact with non-sterile environment.For the above hypothetical patient, we see that with no G-CSF support, or with daily injections, n(t) < n trc for substantial periods.Hence, even the tiniest bacterial contamination can cause acute infection.For the patient receiving bi-daily injections, having nadir neutrophil values of N nadir ≈ 8.7 × 10 5 , we find from equation ( 16) that the critical bacterial influx sustained by this patient is s c ≈ 306#bac/minute whereas numerically we find that s c ≈ 400#bac/minute.
Notably, this critical influx value depends sensitively not only on the neutrophil counts at the nadir, but also on the neutrophils ability to defeat the bacteria, namely the kill term parameters.Indeed, when we increase the kill parameter α by 20% this patient does not develop acute infection under the optimistic model even when the bacterial flux is increased to s = 3500#bac/minute (magenta curve in Figure 4e).The estimated asymptotic critical bacterial influx for this patient that has more efficient neutrophil kill rate is s c ≈ 4250#bac/minute whereas numerically we find that actually s c ≈ 3750#bac/minute.
Figure 5 shows how this critical bacterial influx value depends on each of the parameters of the model ( 13) (when all others are held fixed).Numerically, we define s c as the smallest influx value for which acute infection (b exceeds 10 8 bac/ml) develops during the 21 days of the chemotherapy cycle.We see that the critical influx value depends very sensitively on α, quite sensitively on ρ, N nadir and η and less so on the other parameters β, δ, γ (notice the different scales on the y axis of Figure 5).Notice that the analytic asymptotic expression for s c (equation ( 16), blue dashed lines in Figure 4d) and the numerics provide similar behavior yet do not agree even for small s values.The main reason for this discrepancy is that the numerical definition of s c (which is the medically relevant definition) does not coincide with the bifurcation value that is approximated by equation (16).In this definition we purposely limit the integration time to 21 days whereas close to the bifurcation point this time is not sufficiently long to get to the 10 8 threshold (this discrepancy indeed persists even when n is fixed, and becomes smaller when the integration time is increased or the acute infection threshold level is decreased).
Notably, similar behavior appears when the constant bacterial influx is replaced by a two-day infection episode at the nadir.On the other hand, infection episode at the onset of the treatment or after the neutrophil recover do not induce infections in this simplified model.
A few remarks regarding this model are in order.First, we should remark regarding the choice of the parameters in Table 1.These are the parameters presented in [15] supplement (Tables 1 and mean values of Table 2), with the following three modifications that were employed to obtain nontrivial results as explained next.First, we set the control parameter, AM C, to 0.2.This choice was made so that with the bi-daily G-CSF treatments, the neutrophils level do exceed the S. epidermidis bacteria critical value n trc = 7.99 × 10 5 .A smaller AMC value leads to the trivial result that no G-CSF treatment may help to avoid the development of acute infection.Notably, this value is larger than the standard medical threshold level N c .Second, the starting day of the G-CSF therapy are slightly altered from those appearing in [15].Instead of starting the G-CSF injections 6 days after the neutrophils decline we begin the treatment 2 days after the first decline.Third, the neutropenia duration is shorter by one day from that appearing in [15].These two last changes are made so that with the 10 days bi-daily injections the patient neutrophil levels do not drop below n trc = 7.99 × 10 5 .Other choices (of a later start of the injections or of more sever or of longer neutropenia with the same AMC value and without additional injections) always lead to the development of an acute infection for very small influx parameter even with the bi-daily injections.
As mentioned before, this combined model neglects many effects.Most of the neglected effects are negative, in the sense that the actual bacterial killing is expected to be lower than the killing accounted by in the model, hence this model can be viewed as optimistic.We list next some examples for such negative effects.The cumulative effect of the phagocytosis leads to an increase in the tissue neutrophils death rate; The neutrophils rate of migration to the infected sites is not immediate; In-vivo chemotaxis may slow down the phagocytosis; The damage caused by the neutrophils influences the function of the incoming ones; Anti-inflammation factors cause a decrease of the neutrophil recruitment into the tissue.On the other hand, three positive effects are neglected as well.The existence of other phagocytes, killer cells and the complement system in the tissues; The natural longer life-span the neutrophils have in the tissue [52]; The recruitment of the adaptive immune response after about a week from the onset of an infection.The two first positive effects may be especially important in the first days of the neutropenia and may help to overcome minor bacterial influx: these other defence mechanisms and older tissue neutrophils may remain on guard even when the blood neutrophils are already at very low levels.The last positive effect may be relevant for slowly evolving long lasting infections of rare bacteria strains and less so for the common bacterial infections of neutropenic patients [53,54].6. Discussion.We first discuss the modeling issues that arise from the mathematical formulation of the one dimensional model of the bacteria-neutrophil interaction (sections 2-4).We then discuss some of the possible medical implications of the new model that is introduced here, the neutropenic in-vivo model (section 5).
The axiomatic approach for modelling the bacteria-phagocyte interaction results in a precise mathematical classification of all possible bifurcation diagrams that may emerge from the very general class of models that describe this interaction.The different types of bifurcation diagrams (in particular the simplest monostable vs bistable diagrams) correspond to different biological behaviors that can be tested in in-vitro experiments.Notably, these behaviors have not been detected prior to the derivation of this simple class of models, yet, retrospectively, as pointed out in [ 6], the bistable behavior does appear in the previously published data [14].Another prediction of the model is that some bacteria strains should also exhibit a monostable regime (e.g. when the nutrient supply is decreased below some threshold).This behavior may correspond to the non-pathogenic bacteria strains that reside in our bodies.While we believe both types of behaviours should indeed exist, the bi-stable one is probably more common in bacteria strains studied in the medical context, as they can produce infections more easily [6].
Our systematic study of the bacteria-phagocyte interaction lead to new observations regarding two models that are often employed as a building block in population-dynamics: models with quadratic kill terms and models with ratiodependent kill terms.In Section 4.2 we showed that one dimensional models with a reasonable natural growth term and a quadratic mass-action kill term cannot admit bistability.This omission occurs even though a formal linearization of the valid model of Section 4.1 at small b produces a model with a mass-action kill term.Mathematically, this is clear.Linearization in b of the kill term amounts to neglecting nonlinear terms of the form nb 2 .Away from the bifurcation point n trc , at small concentrations, one can certainly neglect such terms.However, near the bifurcation, the magnitude of these terms determines which type of bifurcation occurs.In Section 4.3 we showed that models that use ratio-dependent kill term have a peculiar behavior at small bacterial-concentration values.In particular, these models are incompatible with the well established biological notion of having a finite positive critical neutrophil concentration that is sufficient to defeat small bacterial infection.Notably, the ratio-dependent point of view was employed by several experimentalists [32,34,33,35].In the ecological modeling community, their is a long standing controversy when these modes are appropriate, see e.g.[55,56,57].Summarizing, these two families of models do not adequately describe the behavior near n trc .For neutropenic patients, the behavior near this bifurcation point is, from a medical point of view, the most important regime.
The implications of this observation regarding higher dimensional models that use mass-action kinetics should be carefully examined.In these multi-dimensional models, the above observations regarding the bifurcation diagrams of the 1d system simply provide the form of the b-nullcline.If the multi-dimensional system has some time separation so that the b dynamics is much faster than the slow phagocytes (or more generally the predator) dynamics 11 , this nullcline serves as a backbone for inducing the slow-fast dynamics.Models that employ quadratic mass-action kinetics produce a monotonically decreasing nullcline that corresponds to a stable manifold of the larger system -namely -in such models the dynamics in b hardly plays a role.On the other hand general models with time separation that admit bistable behavior are expected to have much richer multi-dimensional dynamics.
As a first step for studying this complex in-vivo multi-dimensional dynamics, we introduced a new intermediate model which we call the optimistic model (section 5).This model provides a simplified view on how a neutropenic patient fights a constant invasion of S. epidermidis bacteria.Here, the neutrophil levels of the parameterised model ( 13) are taken to be the neutrophil levels in the blood of a chemotherapy patient, as described by equation (19) (following [15]).We note again that many effects are neglected by such a naive model.We expect that these effects will mainly reduce the neutrophils ability to fight the bacteria.Namely, we believe that typically the neutrophils blood level exceeds the tissue level and the kill parameters (α, 1/γ, 1/η) slowly decrease with the time of the infection.Thus, when this combined model predicts the development of an acute infection within a couple of days, we expect that the patient will be indeed in a great danger of developing such an infection in 4-5 days.On the other hand, when this simplified over-optimistic model predicts slower development of an infection or a recovery we cannot rely on it with regards to the true in-vivo dynamics.
This model demonstrates that the bistability of the bacteria-neutrophil dynamics leads to large variability in the response of neutropenic patients to the degradation of their barriers.It also demonstrates that the medical notion of a threshold neutrophil level is sensible, though the universal value of N c needs to be individualized.For a given parameter set and a given bacterial influx level, this threshold corresponds to the saddle node bifurcation value n sn (s).When the neutrophils plummet below this value for too long acute infection develops.In particular, if n drops below n trc even the tiniest bacterial influx causes an infection.Notably, for the S. epidermidis bacteria, this critical value is larger than N c , the medically accepted averaged threshold.This finding may perhaps be reconciled by introducing additional tissue kill mechanisms into the model.
Put differently, we have seen that for a given neutropenic profile with a minimal neutrophil value n = N nadir , there is a critical bacterial influx s c (n) so that beyond this value the patient develops acute infection.We demonstrated that this threshold bacterial influx depends most sensitively on the kill term parameter α and also quite sensitively on ρ, N nadir and η, whereas its dependence on the other parameters is rather mild.While the dependence on ρ and N nadir is not surprising and is indeed in line with current treatment protocols (ρ can be decreased by antibiotics and N nadir is increased by G-CSF injections) the extremely sensitive dependence on α is a new finding that may have new implications on diagnostics and treatments.We believe that further exploring the value of s c in this model and in its extensions may provide important medical information.Potentially, after more clinical data is gathered regarding the common values of s and of the kill term parameters, it can tell us whether, for a given treatment protocol, a patient will be able to overcome the body's internal bacterial population and whether the patient can withstand contact with non-sterile environment.Lemma 6.1.For s = 0, the origin is a fixed point for all n, and there exists n trc > 0 such that b = 0 is unstable for n < n trc and stable for n > n trc .
By A1 F is C r , r > 2, and we showed that ∂ b F (0, 0) > 0, for sufficiently small n ∂ b F (0, n) > 0. However, by equation ( 8) of assumption A5, for small b > 0 and for ñ > n 0 , F (b,ñ)−F (0,ñ) b < 0 and hence ∂ b F (0, ñ) ≤ 0.Moreover, the monotonicity of the local kill rate in n (equation (7) of assumption A4) implies that for all n > ñ we have ∂ b F (0, n) < 0. Therefore, by the Intermediate Value Theorem, there exists a point (n trc ), such that ∂ b F (0, n trc ) = 0.Moreover, by the monotonicity condition of Eq. ( 7) this bifurcation point is indeed unique, dividing the positive n axis to unstable and stable regimes.Lemma 6.2.At n = 0, at sufficiently small s, equation (1) has a unique positive equilibrium solution and it is stable.
Proof.: From the limited growth model (assumption A3), and in particular from equation ( 5), it follows directly that at n = 0 there is a unique positive fixed point, b = b f (s).
Moreover, (5) implies that for all > 0, and thus, generically (specifically, using equation (10) of A6), we get that ∂ b F (b f (0), 0, 0) < 0. Therefore, for sufficiently small s this strict inequality persists and b f (s) is stable.This Lemma shows that while there is a single branch of positive solutions emanating from ((b f (s), 0, s) when n is found as a function of b (i.e. this branch (b, n * (b; s), s) is unique for b > 0), there can be several branches of positive solutions (b * (n; s), n, s) that may meet at the bifurcation points.Next we establish that the stability of the fixed points along such a branch and its monotonicity property are directly related.
Notice that in models with mass action kill term (see section 4.2) ∂ 2 ∂b 2 F (0, n trc , 0) = ∂ 2 ∂b 2 (B(b, s) − D(b))| (b,s)=(0,0) .The above is negative in commonly used models of the natural bacterial growth like the logistic and the Gompertz [31] models.We conclude that when such natural bacterial growth models (with negative second derivative at the origin) are combined with the mass-action kill term, only the supercritical normal form (minus sign in (21)) can be realized.
Let n * (b; 0) denote the unique solution of equation ( 1) at s = 0 that emanates from (b f , 0, 0) (recall lemma 6.3).We next show that this branch connects to the n axis exactly at the bifurcation point: With the above preparation we are now ready to prove the first theorem: Proof.[of Theorem 3.1] By Lemma 6.6 the only intersection of n * (b) with b = 0 is at n = n trc .If ∂ 2 F ∂b 2 (0; n trc , 0) < 0 by Lemma 6.5, the intersecting branch is stable, thus, starting from the stable branch at n = 0 there must be an even number of generic bifurcation points to get a stable branch intersecting the b = 0 branch.If ∂ 2 F ∂b 2 (0; n trc , 0) > 0, by Lemma 6.5, the intersecting branch is unstable, thus, starting from the stable branch, there must be an odd number of generic bifurcation points to get an unstable branch intersecting the b = 0 branch.
Finally, we note that there cannot be any additional branch of fixed points, namely the equation F (b, n, 0) = 0 cannot have any additional non-trivial (non-zero b) solution curves.By Corollary 2 there cannot be any closed fixed point solution curves.We now show that any other option will lead to a contradiction to one of the observations regarding the solution curves properties.Additional intersections of the solution curves with the n-axis (respectively b-axis) contradict Lemma 6.1 (respectively 6.2).Asymptotic branch to the n axis is impossible since it implies that either a stable branch (by Lemma 6.4) asymptotes a stable branch or that an infinite sequence of bifurcations occur as b → 0, contradicting Lemma 6.5.Finally, asymptotic branch or an infinite sequence of bifurcations that occur as n → 0 are ruled out by A3 and continuity (A1).
Proof.[of Theorem 3.3] We first note that for sufficiently small s, all the zero flux saddle node bifurcation points that occur at positive bacteria concentrations persist (beyond the general persistence arguments for the saddle node bifurcation, we notice that (4) implies that ∂F ∂s (b, n, 0) > 0 at these bifurcation points so for small s these bifurcation points simply shift to the right, see also Lemma 3.2).Thus to complete the proof we only need to examine the behavior near the zero flux transcritical bifurcation point.By lemma 6.5 for s = 0 the model is diffeomorphic to the normal form ẋ = ±x 2 − μx, at (b * , n * ) = (0, n trc ) where the bifurcation point is shifted in the normal form to the origin.Removing the constraint that b = 0 is a solution by introducing a small influx s is equivalent to introducing a symmetry breaking parameter into the normal form: ẋ = ±x 2 − μx + , where, by equation ( 4), > 0.
In particular, case (1) of Theorem 3.1 corresponds, for sufficiently small s and for (b, n) values near (0, n trc ), to the transcritical unfolded normal form ẋ = −x 2 − μx + , > 0. This unfolded form has no bifurcation point near (x, μ) = (0, 0).Instead, it has a monotone decreasing stable positive branch, namely, the transcritical bifurcation dissolves into a regular point.Thus, by Theorem 3.1(1) and the above remark we indeed have here an even number (possibly zero) of saddle node bifurcations.
Similarly, in case (2) of Theorem 3.1, for small positive s the model is diffeomorphic to the unfolded normal form ẋ = x 2 − μx + , > 0. This normal form has two saddle-node bifurcations near (x, μ) = (0, 0), yet only one of these corresponds to positive x (and thus b) values.Thus, the transcritical bifurcation point transforms to a single additional saddle node bifurcation point.Therefore, by Theorem 3.1(2) and the above remark, for sufficiently small s there is again an even number of saddle-node bifurcations and no other types of bifurcations may emerge.In particular, by corollary 1, at s = 0 we have here at least one saddle-node bifurcation at positive b value.Thus, for small s, there are at least two saddle-node bifurcation points.

A 1 .
The functions B(b; s), D(b), D K (b; n) are sufficiently smooth (C r , r > 2) and all the birth and death rates 5 are bounded.

d 2 F 2 > 0 Figure 1 .
Figure 1.The bifurcation diagram presents the bacteria equilibrium dependence on the neutrophils.The red region indicate prosperity for initial bacterial population, and blue region indicate a decrease in the initial bacterial population.Solid lines indicate a branch of stable fixed points, and dashed lines indicate a branch of unstable ones.(a)Monostable: in this model, if n < n trc , regardless of the initial concentration of the bacteria (other than zero), the bacteria concentration converges to a positive bacteria equilibrium.For n > n trc the bacteria concentration always converges to the origin.(b) Bistable: in this case there are two critical points, n trc and n sn .For n < n trc the bacteria concentration converges to the positive stable branch.For n > n sn the bacteria concentration converges to the origin.When n trc < n < n sn we have bistability, so the final state of b depends on whether the initial bacteria concentration is above or below the unstable fixed point (dashed line).The presented bifurcation curves correspond to the explicit model (13) at (ρ, δ, α, γ, η, s) = (10 −2 , 10 −3 , 10 −8 , 10 −8 , 10 −10 , 0) with β = 9 • 10 −9 (a) and β = 10 −9 (b).monotonic function as claimed.Similarly, by A3-A4, fixing b = b * > 0 and looking for a solution of the form F (b * ; n * (s), s) where F (b * ; n * (0), 0) = 0, we obtain that dn * (s) ds = − ∂sF ∂nF = ∂sB(b;n,s) ∂nD K (b,n) > 0. This inequality holds for all positive equilibria solutions of the zero flux model (b * ; n * (0)), including the bifurcation points.

: d 2 F 2 > 0 Figure 2 .
Figure 2. Bifurcation diagrams with zero and small positive bacterial influx.Bifurcations curves of (13) at s = 0, 2 • 10 5 , 4 • 10 5 are shown in black, magenta and green, respectively (all other parameters are as in Fig 1).The bifurcation curves are shifted to the right as s is increased.(a) Monostable: the zero s transcritical bifurcation point at n trc disappears when s > 0. (b) Bistable: the zero s transcritical bifurcation becomes a saddle-node bifurcation that appears at a distance O( √ s) from the transcritical bifurcation point.(Insert) A close-up around the transcritical bifurcation of (b), which shows that when s > 0, the transcritical bifurcation becomes a saddle-node bifurcation.
Increased s and α

Figure 4 .
Figure 4. Bacterial dynamics during a neutropenic cycle according to equation(19).The dynamics of (a) bacteria, (b) neutrophils, and (c) G-CSF for patients with bacterial influx of s = 100 [#bac/ml/min], undergoing chemotherapy treatment with 3 different G-CSF support protocols.No treatment -cyan curves, daily injections -blue curves and bi-daily injections -green curves.In (d) these trajectories are projected onto the neutrophils-bacteria plane.In (e) the response of the patient receiving bi-daily injections to an increased bacterial influx is demonstrated.When the influx is increased from s = 100 (green, as in (a)) to s = 400 (red) acute infection develops.When the kill term parameter α is increased by 20%, a bacterial influx as large as s = 3500 [#bac/ml/ min] results in an infection that is under control (magenta curve).

Figure 5 .
Figure 5. Sensitivity of the critical bacteria influx of a patient receiving bi-daily G-CSF injections to parameter variations (neutrophil counts are as in the green curve of Fig 4a).Red curve: numerical values of the critical flux by day 21, Blue curves: the asymptotic expression for s c (Eq. (16)).Black points: the parameter values used in Fig 4.

Lemma 6 . 3 .
For any fixed s ∈ [0, s max ), the set of positive fixed points of (1), {(b, n, s)|F (b, n, s) = 0, b > 0}, contains the connected set {(b, n, s)|n = n * (b; s), b > 0} where n * (b; s) is the unique solution of the equation F (b, n * (b; s), s) = 0 that passes through the positive b axis (satisfying n * (b f (s); s) = 0).Moreover, the only points at which the function n = n * (b; s) is not locally invertible with respect to b are the bifurcation points of (1).Proof.: By A4, for all b > 0, ∂ n F (b, n, s) = ∂ n D K (b, n) > 0, and thus near any fixed point of (1) n * (b; s) may be uniquely found by the Implicit Function Theorem.In particular, since by A3 for any fixed s there is a unique point b f (s) such that F (b f (s), 0, s) = 0, we obtain that indeed the intersection of the curve n * (b; s) with the positive b axis must occur at b f (s) and that the branch emanating from n * (b f (s); s) = 0 is uniquely defined and finite for all b ∈ (0, b f (s)] (recall that ∂ b F (b, n, s) is finite in the domain Ω).Similarly, the function b * (n; s) that solves F (b * (n; s), n, s) = 0 is defined by the implicit function theorem near any fixed point (b, n, s) s.t.∂ b F (b, n, s) = 0 (i.e., at all fixed points that are not bifurcation points).Moreover, at such points, by the inverse function theorem, the inverse with respect to b of the function n = n * (b; s) is simply b * (n; s).

Lemma 6 . 4 .
If b * (n; s) > 0 is a branch of fixed point of (1) which maintains its linear stability for all n ∈ [n − , n + ], then it is monotone in n; If b * (n; s) is stable in this interval then it is decreasing in n and if b * (n; s) is unstable then it is increasing in n.Proof.By assumption F (b * (n; s), n, s) = 0 and ∂ b F (b, n, s)| b * (n;s),n,s = 0 for all n ∈ [n − , n + ], therefore, by implicit differentiation ∂b * (n;s) ∂n

Figure 6 .Corollary 3 .Lemma 6 . 5 .
Figure 6.The model assumption exclude bifurcations leading to an increasing stable branch or decreasing unstable branch.

Table 1 .
(see Table1for the parameters, dimensions and values), here The chosen Neutrophils-G-CSF model parameters.The duration of the chemotherapy, the AMC value and the starting day of the G-CSF therapy are slightly altered from the n,s).Thus, by A4 and since ∂ b F (b * (n; s), n, s) has a fixed sign here, b * (n; s) is indeed monotone on this segment.Moreover, since the numerator is positive (by A4) it is decreasing in n for stable branches (where ∂ b F (b * (n; s), n, s) < 0) and increasing for unstable branches (∂ b F (b * (n; s), n, s) > 0).