A computational model to inform effective control interventions against Yersinia enterocolitica coinfection

: The complex interplay among pathogens, host factors, and the integrity and composition 1 of the endogenous microbiome determine the course and outcome of gastrointestinal infections. 2 The model organism Yersinia entercolitica (Ye) is one of the ﬁve top frequent causes of bacterial 3 gastroenteritis based on the Epidemiological Bulletin of the Robert Koch Institute (RKI) published on September 10, 2020. A fundamental challenge in predicting the course of an infection is to understand 5 whether co-infection with two Yersinia strains differing only in their capacity to resist killing by 6 the host immune system may decrease the overall virulence by competitive exclusion or increase 7 it by acting cooperatively. Herein, we study the primary interactions among Ye, the host immune 8 system and the microbiota, and their inﬂuence on Yersinia population dynamics. The employed 9 model considers two host compartments, the intestinal mucosa and lumen, commensal bacteria, 10 the co-existence of wild-type and mutant Yersinia strains, as well the host immune responses. We 11 determine four possible equilibria: the disease-free, wild-type-free, mutant-free, and co-existence of 12 wild-type and mutant equilibrium. We also calculate the reproduction number for each strain as a 13 threshold parameter to determine if the population may either be eradicated or persist within the 14 host. We conclude that the infection should disappear if the reproduction numbers for each strain fall 15 below one, and the commensal bacteria’s growth rate exceeds the pathogens’ growth rates. These 16 ﬁndings will help inform medical control strategies. The supplement includes MATLAB source script, 17 Maple workbook, and ﬁgures. 18

. The variable symbols and their meaning. The lumen is abbreviated with L, the mucosa as M.
We also indicate each variable with mut or wild-type (wt) to denote to which population they refer. An upper-case I refers to the immune system. Upon oral co-infection, a Ye population enters the small intestinal lumen. Due to their particular 83 virulence traits, some Ye cells are also able to enter the mucosal site.

84
As shown in the model eqs.
(1) to (7), the population dynamics of bacterial species consist of 85 both growth, i.e., an increase of populations due to a distinct growth rate α, and a reduction through 86 peristalsis, which moves the bacterial populations towards the colon where they finally end up in feces.

87
Peristalsis will only influence the populations in the lumen. Furthermore, bacterial populations in the 88 mucosa can additionally be reduced through killing by the host immune attack.

89
Both sites, the lumen and the mucosal site, are colonized by commensal bacterial species but to a 90 different extent. The colonization capacity of the lumen considerably exceeds that of the mucosal site.

91
This is due to the host's natural mechanisms to limit bacterial contact to the epithelium through physical 92 barriers, like the adherent mucus layer and a high concentration of anti-microbial peptides (AMPs).

93
Hence, we assume a much lower carrying capacity of the mucosal site compared to the luminal 94 compartment. Another assumption is that as soon as bacterial numbers exceed the mucosa's carrying 95 capacity, they spill-over to the luminal site, feeding the luminal populations [10]. 96 The model's immune system represented by I (see table 1) is activated as soon as at least one  In contrast to commensal bacteria, Ye exhibits a number of virulence traits to evade the host immune response. This is implemented in our model by different immunity adjustment factors for either the wild-type strain or the different mutant strains. Consequently, the number of commensal species is more affected by the host immune attack than the number of Ye mutant strains, which are, of course, more affected than the wild-type strain. Our model's final output is the number of bacteria or colony-forming units (CFUs), finally ending up in feces. These population dynamics are represented as the following a system of differential equations [10]: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 October 2020 doi:10.20944/preprints202010.0386.v1 For a detailed description of the parameters, see table 2. Besides the seven differential equations describing the population dynamics of the Ye strains, the commensal bacteria and the immune system at the mucosal site, and all three populations (commensal bacteria, wild-type Yersinia, and mutant Yersinia) at the luminal site, the model consists of different mass-action terms describing the spill-over rates from the mucosa into the lumen and describing the growth rates of the bacterial populations, which are limited through a given capacity "C" and are defined as follows: and model parameters are shown in table 2. Since the system is the modeling of Yersinia population in The basic reproduction ratio R 0 is calculated by the fraction of the transmission rate for each 105 strain (spill-over) and the average infectious period for each strain in compartments as defined by In addition, van den Driessche and Watmough [25] defined R 0 as a general compartmental disease transmission model suited to heterogeneous populations that can be modeled by a system of ordinary differential equations. The authors, thus, divided systemẊ = f (X) as two different main compartments with the first compartment F i (x) defined as the rate of the appearance of new infections in compartment i, V + i as the rate of transfer of individuals into compartment i by all other means, and V − i as the rate of transfer of individuals out of the compartment i where can be written . . , n. Therefore, the disease transmission model consisting of non-negative initial conditions will be as follows: It is assumed that functions F i (x) and V i = V − i + V − i are continuously differentiable at least 108 twice in each variable, and they satisfy the assumptions items A(1) to A(5) described below [25,26]: A(5) The systemẏ = G(0, y) has unique asymptotically stable equilibrium, y * 114 Theorem 1. Consider the disease transmission model given byẊ = f (X) with f (X) satisfying conditions 115 items A(1) to A(5). If X 0 is a disease-free equilibrium (DFE) of the model, then X 0 is locally asymptotically 116 stable if R 0 < 1, but unstable if R 0 > 1, where R 0 acts as a threshold parameter.

118
Several models found in the literature have been used to show that when R 0 crosses the threshold, R 0 = 1, it can act as a bifurcation parameter and transcritical bifurcation takes place. That is, local asymptotical stability is transferred from the disease-free state to the new (positive) equilibria. In ordinary differential equations, we will encounter bifurcations of equilibrium and periodic orbits, which are typical in the sense that they occur when a small smooth change is made to the threshold parameter values (the bifurcation parameters µ) of a system. It causes a sudden qualitative or topological change in the behavior of the system. Consider the continuous dynamical system described by the Ordinary Differential Equation (ODE) A local bifurcation occurs at (X 0 , µ 0 ) if an eigenvalue with zero real part is included in the Jacobian 119 matrix of the system. If the eigenvalue is equal to zero, the bifurcation is a steady-state bifurcation.

120
Therefore, we now recall the analysis of the center manifold near the critically (X = X 0 , R 0 = 1), 121 which allows clarifying the direction of the bifurcation near the bifurcation point using the Center 122 Manifold Theorem 2. This theory describes the determination of the local stability at the non-hyperbolic 123 equilibrium (linearization matrix has at least one eigenvalue with zero real part) and the existence of 124 another equilibrium (bifurcated from the non-hyperbolic equilibrium).

126
The center manifold theorem provides a means for systematically reducing the dimension of the 127 state spaces which need to be considered when analyzing bifurcations of a given type.

128
Theorem 2 (Center Manifold Theorem for Flows). Let f be a C r vector field on R n vanishing at the origin f (0) = 0 and let A = D f (0). Divide the spectrum of A into three parts, σ s , σ c , σ u with Let the (generalized) eigenspaces of σ s , σ c and σ u be E s , E c and E u , respectively.

133
To achieve this, consider the general systeṁ where Jx is the linear part of the system. We must first find the differential equations on its center 134 manifold and then reduce the system to its normal form. Without loss of generality, we assume that 135 x = 0 is the fixed point of interest for the system.

136
Suppose J has n c eigenvalues with zero real-part and n s eigenvalues with negative and positive real-part and n = n c + n s . Using the eigenvectors of J to form a transformation matrix, the system can be rewritten in block matrix form as where (x c , x s ) ∈ R n c × R n s , A ∈ R n c × R n c and B ∈ R n s × R n s . With the eigenvalues of zero real-part, the Center Manifold Theorem 2 guarantees that there exists a smooth manifold W c = {(x c , x s ) | x s = q(x c )} the equilibrium point such that the local behavior in the center direction of the system is qualitatively the same as that on the manifold. By differentiating x s = q(x c ) , we geṫ x s = Dq(x c )ẋ c . Substituting eq. (19) into the previous identity and rearranging the equation, we get By solving for q(x c ), we get a function describing the center manifold. In general, q(x c ) cannot be solved explicitly. Instead, substituting a Taylor expansion q(x c ) = ax 2 c + O x 3 c into eq. (20), we can find the coefficients for the expansion by balancing the lower order terms. Based on q(x c ), we now have a system in the reduced form:ẋ The following proposition helps us understand the type of transcritical bifurcation is forward or

147
We used the model in eqs.
(1) to (7) to understand how the dynamics change following variation 148 of different parameters. We calculated the equilibria of eqs.
(1) to (7), conduct a linear stability analysis, 149 and identified the analytical conditions that lead to a transcritical bifurcation. For mathematical convenience, we divide the model in eqs.
(1) to (7) such that the first four equations correspond to infected individuals. The distinction between infected and uninfected populations must be determined from the model's epidemiological interpretation and cannot be deduced from the structure of the equations alone. Thus, we have while we assume that the growth rate α (B) of the endogenous commensal bacteria is higher than the 152 Ye growth rates α (wt) and α (mut) , respectively.

153
On the other side, the Yersinia model has three compartments (mucosa, lumen, immune system), which are analyzed separately. The model system is analyzed in a suitable, feasible region. All forward solutions of the system are feasible ∀t ≥ 0 if they enter the invariant region The existence of equilibrium points for system eqs.

154
• The trivial equilibrium point is as an origin equilibrium (0, 0, 0, 0, 0, 0, 0). This solution appears 155 when all populations are extinct. For all parameters, this point never becomes stable due to 156 positivity of eigenvalues in A2.
• The first equilibrium point appears in the absence of Yersinia Y (1) to (7) has a disease-free equilibrium, which is given by and It describes a disease-free state whereby only the commensal bacteria persist. In order for the 158 disease-free state P 0 to be biologically meaningful, the conditions γ < α (B) and β < α (B) must 159 hold. These conditions correspond to the maximal growth rate of intestinal bacteria exceeding 160 the rate at which intestines are charged and the maximal immunity action, which is not that 161 strong in the absence of Yersinia strains. However, the population of the immune system is in a 162 maximum of its carrying capacity (in health, not in fighting with any infection).

163
• A second equilibrium corresponds to the commensal bacteria's persistence and the Yersinia mutant strain in the absence of the wild-type strain. Without loss of generality, the commensal bacteria are supposed to be zero because they are not infective. This point is obtained by setting • The other equilibrium corresponds to the persistence of commensal bacteria and the Yersinia wild-type strain in the absence of the mutant strain. Without loss of generality, the commensal bacteria are supposed to be zero because they are not infective. This point is obtained by setting Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 October 2020 doi:10.20944/preprints202010.0386.v1 • Finally, the last equilibrium point corresponds to a state of the co-existence of wild-type and mutant Yersinia strains. This point is achieved by supposing B M = B L = 0: to make the equilibrium point biologically 164 meaningful. (1) to (7) is given by the equation A1 in the appendix.

169
The disease-free equilibrium of the model is the steady-state solution of the model in the absence of the disease. The eigenvalues of the Jacobian matrix A1 at this point were calculated as follows: in which it is assumed that α (wt) < α (B) and α (mut) < α (B) . Based on the general compartmental model 170 describing an infectious disease transmission within a heterogeneous population, the host population 171 is grouped in two general classes, the infected and uninfected compartments. Therefore, the system 172 eqs.

173
In classical epidemic models, it is common to observe that a disease-free equilibrium loses its 174 stability for R 0 = 1, and a transcritical bifurcation occurs. A transcritical bifurcation is when a fixed can appear. We mathematically analyze this aspect within the structure of the model to assess which 181 parts of the structure might be responsible for the occurrence of the transcritical bifurcation.
182 Let us consider system eqs.
(1) to (7) with the above-calculated eigenvalues. Using the eigenvectors A10 shown in the appendix as a basis for a new coordinate system, we set the transformation matrix whose columns are the eigenvectors; Under this transformation, the system eqs.
(1) to (7) so that the linear part is now in a standard (diagonal) form. In the (u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 ) coordinates, the center manifold is a curve tangent to the u 1 − axis. The projection of the system onto the u 1 − axis is obtained by setting u 2 = u 3 = u 4 = u 5 = u 6 = u 7 = 0 in the equation foru 1 . It yieldsu 1 = 0.
Since E c is one dimension, we can approximate h i (u 1 ) by Taylor expansion such that u i = h i (u 1 ), i ∈ {2, . . . , 7}, satisfy in the following equationṡ Thus, the reduced system isu The advantages of a center manifold are clear from this calculation. We may study a one-dimensional  To understand the role of basic reproduction number, we also define a single reproduction number for commensal bacteria R com 0 computed as an expected number of a secondary case of commensal bacteria produced by a single commensal bacteria: Due to biologically meaningful of disease-free state P 0 and holding the conditions β < α (B) and for R wt 0 < 1 and R mut 0 < 1, and the uninfected state is stable (not asymptotically stable) as none of the strains persist. As soon as the two equilibria collide non-destructively, exchanging their stability and 208 resulting in the Jacobian matrix having a single eigenvalue that is equal to zero, then either the mutant 209 strain or the wild-type appears and persists for R mut 0 or R wt 0 , respectively.

210
As shown in fig. 1, when the disease-free equilibrium loses its stability, different scenarios can 211 occur. Herein, we discuss the different cases: Bogdanov-Takens bifurcation [29,30]. The disease-free equilibrium P 0 is loosing its stability and 217 the wild-type-free and mutant-free are including one simple zero eigenvalue (λ 3 = 0), meaning 218 that the dynamic of the model are changing as the target parameter is in the threshold value.

219
(II) If R wt 0 > 1, the wild-type strain equilibrium in the region II will persist when R mut 0 < 1. The 220 wild-type strain will spread and possibly persist within the host population. In general, for a 221 strain to persist, its basic reproduction number has to be strictly greater than one. Therefore, in 222 this region, the disease-free, mutant strain and co-existence state exchange stability: P 0 becomes 223 unstable, P 2 becomes locally asymptotically stable, and P 1 and P 3 remain unstable. This means 224 that the immune system could kill one of the strains more efficiently.

225
(III) If R mut 0 > 1, the mutant strain equilibrium in the region III will persist when R wt 0 < 1. The mutant 226 strain will spread and possibly persist within the host population since its basic reproduction 227 number is greater than one. Therefore, in this region, the disease-free, wild-type strain, and 228 co-existence state exchange stability: P 0 becomes unstable, P 1 becomes locally asymptotically 229 stable, and P 2 and P 3 remain unstable. This means that the immune system could defeat the 230 wild-type strains. However, the risk of happening this situation is low because the mutant! 231 strains are influenced more efficiently than wild-type strains by immune action.

262
Following the analysis of the model in respect to α (B) , two cases are considered. The first case is 263 when α (wt) and α (mut) are equal. In this case, we do not have co-existence of wild-type and mutant (the and Y 3(mut) L will be equal to zero). The second one is when α (wt) and α (mut) have 265 different values as the co-existence equilibrium is achieved and be biologically meaningful.

269
Through the definition of R com 0 , R wt 0 and R mut 0 , we display a x-scale as α (B) fig. 3. This shows the 270 appearance of different equilibrium of model eq. (14) in terms of the basic reproduction numbers.

271
As shown in fig. 3, we expect no co-existence equilibrium when α (wt) = α (mut) and only wild-type 272 equilibrium exists (mutant-free equilibrium) when 1.76 < α (B) < 11.20. However, the disease-free The experimental data obtained in the lab has been shown that the Ye mutant strain is eliminated 291 more efficiently than Ye wt at a mucosal compartment. However, since GF mice lack a microbiome, 292 both the Ye wt and the mutant strain can colonize the gastrointestinal tract (GIT) at high numbers. In 293 the SPF-colonized Myd88! (Myd88!) mice, the immune response is weaker [10]. Therefore, as shown 294 in fig. 6, fig. 7, fig. 8, and fig. 9, the following results conclude: fig. 6d and fig. 7d shows a slow reduction of the wt strain in compare with mut strain as κ is 296 changing.

297
• fig. 6b, fig. 7b, fig. 8b, and fig. 9b in compare with similar situations in SPF and MyD88 -/shows 298 elimination of wt and mut strains in less efficient.

299
• fig. 6c, fig. 7c, fig. 8c and fig. 9c shows that wt and the mutant strains increased faster. Besides. 300 the influence κ on wt and the mutant strains can not project properly ( fig. 8f, fig. 9f) in a same 301 speed of producing strains.

302
Our results allow us to state that if mice possess a healthy microbiota and immune system, as 303 long as the growth rate of commensal bacteria is larger than the growth rate of wild-type and mutant 304 Yersinia, then the infection will not spread as Ye strains cannot enter the lumen compartment. Therefore, This region is not biologically meaningful. The only equilibrium in this region is trivial solution.
The appearance of equilibrium point corresponds to a state of co-existence of wild-type and mutant Yersinia strains.
The appearance of equilibrium corresponds to the persistence of the Yersinia wild-type strain in the absence of the mutant strains.
The appearance of free-disease equilibrium point ≈ 2.28 (α (wt) = α (mut) ) of co-existence of wild-type and mutant. The condition of (α (wt) ≠ α (mut) )is necessary for the existence of co-existence equilibrium.
In this case, the region behind the blue point is not biologically meaningful. We do not have a quilibrium point corresponds to a satate > 1 and R wt 0 < 1. Whenever R mut 0 > 1, R wt 0 is already above one; this facilitates the situation 317 for a co-existence equilibrium. This is where both strains persist, exists when both reproductive 318 numbers are above a certain threshold. However, we could not analytically solve stability criteria for 319 the equilibrium due to complexity, but simulations show that this stability criteria exist.  Figure 4. The sensitivity of parameter in respect to α (B) for 336 hours. 0 < α (B) < 1, no biologically meaningful region (out of our interest). α (B) > 1, all populations appear. When 1 < α (B) < 2.28, the region is a region for the appearance of the co-existence equilibrium, but the hypothesis of the co-existence equilibrium is not satisfied. Therefore, wild-type strain does not grow, and the mutant strain is going down slowly. When α (B) > 2.28, this is a region of the appearance of wild-type equilibrium (R wt 0 > 1). Thus, a and e are increasing fast and staying in the maximum level as a and d are going back to zero. Additionally, c and d do not grow anymore (R mut 0 < 1). population-level dynamics, but also suggest that this would require strong non-monotone effects in 333 the immune response to infection.

334
Furthermore, this also raises the question of whether the infection becomes persistent or is 335 resolved. The answer to this question depends on the magnitude of the basic reproductive number, 336 R 0 , since the calculation of the stability of equilibrium is very complicated. As shown, the disease-free 337 equilibrium is locally (asymptotically) stable when 0 < R 0 < 1 and unstable if R 0 > 1. In other 338 words, when 0 < R 0 < 1, every infectious strain will cause less than one secondary infection, and 339 hence the disease will die. When R 0 > 1, every infectious strain will cause more than one secondary 340 infection, and hence Yersinia infection will occur. All public health control measures can usually be 341 based on methods that, if practical, would lower R 0 to below unity. On the other hand, the co-infection 342 equilibrium is locally stable when R 0 > 1 and unstable when 0 < R 0 < 1. Typically, for Ye, we do 343 not observe the spread of infection. However, Ye is a great model system and can be used to predict 344 the spread of infections by more clinically relevant bugs, e.g., enteropathogenic or enterohemorrhagic 345 Escherichia coli. To control the spread of infection, efforts must be made to ensure that the co-infection 346 equilibrium is unstable, i.e., 0 < R 0 < 1.

347
In a co-infection setup, we found two thresholds for the model, R wt 0 and R mut 0 . The threshold 348 associated with the wild-type and mutant strains determine the existence and the stability of the 349 equilibrium point, and we considered four equilibrium points for the model.

350
As the theoretical work agrees with the experimental work, our results may be used to control 351 the spread of infection by other pathogens of the gastrointestinal tract, e.g., enteropathogenic or (a) Co-existence of wild-type and mutant Yersinia strains in the mucosa (b) Co-existence of wild-type and mutant Yersinia strains in the lumen Figure 5. The diagram displays the appearance of the co-existence of wild-type and mutant Yersinia strains as α (B) is changing in a region where R wt 0 > 1 and R mut 0 > 1. a An immune reaction influences the wild-type and mutant Yersinia strains in the mucosa. Additionally, another part of them is spilled over and move to the lumen compartment. Therefore, they are simultaneously increasing to reach their maximum values and decreasing. b The wild-type and mutant Yersinia strains in the lumen are simultaneously increasing to reach the carrying capacity of the lumen compartment. In both a and b, when α (B) is reaching 4.90, the co-existence of wild-type and mutant Yersinia strains is disappearing. enterohemorrhagic E. coli. Furthermore, an additional factor to take into account is whether antibiotics Appendix A.2. The eigenvalues of trivial equilibrium point 372 The eigenvalues for trivial equilibrium point are as follows: Appendix A.6. The basic reproduction numbers 376 The calculation of the basic reproduction number is as follows: The projection of disease-free state Figure 10. Diagram displaying the disease-free state by fixing the parameter values Ye SPF wt/A0 table 3 in model eq. (14). As long as the α (B) is large enough, the disease-free state persists. However, by the reduction in α (B) , the basic reproduction number R wt 0 reaches its threshold value (R wt 0 = 1). This causes changes in the dynamic behavior of the model eq. (14) as shown in fig. 3.