Computational Model Informs Effective Control Interventions against Y. enterocolitica Co-Infection

Simple Summary Medical control strategies for infectious diseases remain enormously important. One germ that can cause gastrointestinal infections is Yersinia enterocolitica. This study investigates and analyzes a computational model to identify the occurrence of disease-free and co-infection states. Thereby, the reproduction number R0 informs us about the germ’s ability to spread disease. Suppose this fundamental quantity takes a value between zero and one. In that case, every infectious strain will cause less than one secondary infection, so the strain will disappear. In contrast, if R0 exceeds one, every infectious strain causes more than one secondary infection, and Yersinia infection strains will persist. A disease-free state occurs when the commensal bacteria’s growth rate exceeds the maximum immune action and the rate at which the intestines release the bacteria. With a large enough commensal bacteria growth rate, this state can be stable. Co-infection occurs when the maximum growth rates of the wild-type and mutant strains become unequal. Studying the immune system’s behavior can result in an infection’s disappearance from hosts with a healthy microbiota immune system. In this case, Yersinia strains do not spread in the lumen when the commensal bacteria’s growth rate exceeds the growth rate of wild-type and mutant Yersinia. Abstract The complex interplay between pathogens, host factors, and the integrity and composition of the endogenous microbiome determine the course and outcome of gastrointestinal infections. The model organism Yersinia entercolitica (Ye) is one of the five top frequent causes of bacterial gastroenteritis based on the Epidemiological Bulletin of the Robert Koch Institute (RKI), 10 September 2020. A fundamental challenge in predicting the course of an infection is to understand whether co-infection with two Yersinia strains, differing only in their capacity to resist killing by the host immune system, may decrease the overall virulence by competitive exclusion or increase it by acting cooperatively. Herein, we study the primary interactions among Ye, the host immune system and the microbiota, and their influence on Yersinia population dynamics. The employed model considers commensal bacterial in two host compartments (the intestinal mucosa the and lumen), the co-existence of wt and mut Yersinia strains, and the host immune responses. We determine four possible equilibria: disease-free, wt-free, mut-free, and co-existence of wt and mut in equilibrium. We also calculate the reproduction number for each strain as a threshold parameter to determine if the population may be eradicated or persist within the host. We conclude that the infection should disappear if the reproduction numbers for each strain fall below one, and the commensal bacteria growth rate exceeds the pathogen’s growth rate. These findings will help inform medical control strategies. The supplement includes the MATLAB source script, Maple workbook, and figures.


Introduction
Forecasting the evolution of infectious diseases is the primary motivation behind the use of mathematical models in biology. Determining those threshold values that predict whether the disease will spread within the host or can be contained is crucial to inform medical control strategies.
Yersinia entercolitica (Ye) is a Gram-negative enteropathogen causing foodborne gastrointestinal infections. Within the small intestine (SI), Ye can adhere to and invade the intestinal epithelial lining mainly via the adhesins Yersinia adhesin A (YadA) [1] and Invasin [2][3][4]. YadA is among the essential virulence factors as a YadA-deficient strain is impaired 97 colonizations, and systemic spread in a mouse model of infection has been shown Upon attachment, Ye can engage its Type III secretion system (T3SS) to inject effector proteins, the Yersinia outer proteins (Yops), into host cells, thus evading the host immune response and establishing a productive infection [5][6][7]. Indeed, the first line of host defense against invading Ye is a massive infiltration of phagocytotic cells. However, Ye can counteract phagocytic killing via its T3SS [6,7]. Together, both the T3SS and YadA contribute to the efficient colonization of the intestinal tract, where Ye induces an inflammatory response that likely accounts for a reduction in density and complexity of the commensal microbiome [8,9]. Several Ye serotypes have been isolated from animal reservoirs and the human gastrointestinal tract, but only a few of them cause disease in humans. Although causing severe distress, gastrointestinal infections are generally self-contained. Typically, healthy persons will only receive symptomatic treatment aimed at avoiding dehydration. However, in individuals with underlying disease, elderly persons, or newborn children, gastrointestinal infections can cause high morbidity and even fatal outcomes. Especially in such patients, identifying those at high risk of developing fatal diseases could help tailor personalized therapeutic interventions to improve the outcome of infection. To this aim, we recently developed a computational model that can calculate Ye population dynamics during gastrointestinal infection and predict pathogen expansion, gut colonization, and infection course [10].
Previously published models have focused on the virulence evolution during co-infection and superinfection by more than one pathogen [11][12][13][14][15][16][17]. Furthermore, Nurtay et al. [18] investigated theoretical conditions in co-existence strains in their research focused on viral populations of wt and mut [18]. Herein, we apply our computational model to predict the behavior of Ye population dynamics during the co-infection of mice by different Ye mutant (mut) strains, lacking distinct virulence factors, and a Ye wt strain in bacterial gastroenteritis. We perform a bifurcation analysis to show how the dynamics of the system change as multiple parameters are varied [19,20].
The existence of a backward bifurcation, i.e., the co-existence of stable disease-free equilibrium with one or more stable endemic equilibria, has significant consequences on the process of producing medical policies aimed at eradicating or controlling an infectious disease within the host [21]. A fundamental parameter in models displaying a backward bifurcation is the primary reproduction number R 0 . R 0 represents the expected number of new infections caused within the host or between hosts when a pathogen enters into the host. If R 0 > 1, after an initial introduction, the infection spreads within the host/or between hosts, thus creating an infectious disease. If R 0 < 1, small initial introductions are not sufficiently transmissible to spread the infection within the host (or between hosts), and the cells get infected within a host. This is called an endemic disease that will fade out. Thus, many control policies like medication/vaccination have focused on reaching coverage levels sufficient to reduce R 0 below 1 [15,16,22,23]. Therefore, in this analysis, we compute the system's equilibrium points, study their stability and the center manifold, and investigate how the different points would affect the infection rate, virulence, and mutation rate. Next, we calculate the primary reproductive number for scenarios in which multiple strains are introduced, and an R 0 number is calculated for each strain. Finally, we discuss the biological implications of our results.

Model Description
Our model for Ye, as described in Geißert et al. [10], considers three different sites with their individual population dynamics. Figure 1 depicts an overview of the model in the form of a Systems Biology Graphical Notation (SBGN) [24] Process Description (PD) diagram [25] based on an illustration by Geißert et al. [10]. The lumen and the mucosal site of the small intestine are illustrated as two separate compartments by abbreviated notations L and M in Table 1. The mucosa and lumen sites include three different population dynamics: commensal bacteria, wt Yersinia, and mut Yersinia.

Mucosa compartment containing Ye
Immune System  Table 1 explains the notations used in the model and this figure.
The black arrows represent processes with an impact on the population dynamics of the entity pools. Arrows pointing from empty set symbols to pool nodes denote an increase in the population size or a decrease if the process arcs point from entity pools to empty sets. Migration across compartments of the respective populations appears as vertical process arcs. Some of these processes receive stimulating effects from the immune reaction or from the size of the Ye populations within the mucosa, as colored arcs indicate. Reference [10] provides a more detailed description of the model's structure. We also indicate each variable with mut or wt to denote to which population they refer. An upper-case I refers to the immune system. SBML [26,27] defines the units item and dimensionless to indicate that a quantity occurs in a piece number or that its unit originates from the cancellation of other units. Besides the bacterial populations, we have the strength of the host immune response. For the sake of clarity, the complex immune cell population dynamics, which is made up of innate and adaptive immunity, are implemented in our model as a single abstract immune response in Table 1. An advantage of this rather abstract immune response is that it is easily adjustable. Only the Ye populations in the mucosa activate the immune response, whereas tolerance to commensal bacteria exists. However, this immune response affects all bacterial populations at the mucosal site. Upon oral co-infection, a Ye population enters the small intestinal lumen. Due to their particular virulence traits, some Ye cells are also able to enter the mucosal site.

Variable Symbol Meaning Units
As shown in Models (1)- (7), the population dynamics of bacterial species consist of both growth, i.e., an increase in populations due to a distinct growth rate α, and a reduction through peristalsis, which moves the bacterial populations towards the colon where they finally end up in feces. Peristalsis will only influence the populations in the lumen. Furthermore, bacterial populations in the mucosa can additionally be reduced through killing by the host immune attack.
Both compartments (the lumen and the mucosal) are colonized by commensal bacterial species but to different extents. The colonization capacity of the lumen considerably exceeds that of the mucosal site. This is due to the host's natural mechanisms to limit bacterial contact to the epithelium through physical barriers, such as the adherent mucus layer, and a high concentration of anti-microbial peptides (AMPs). Hence, we assume a much lower carrying capacity of the mucosal site compared to the luminal compartment. Another assumption is that as soon as bacterial numbers exceed the mucosa's carrying capacity, they spill-over to the luminal site, feeding the luminal populations [10].
The model's immune system represented by I (see Table 1) is activated as soon as at least one Ye cell enters the small intestinal mucosal compartment. This immune activation increases with the number of Ye cells in the mucosa. By triggering the killing of all bacterial cells at the mucosal site, this leads to a reduction in bacterial populations spilling over to the lumen [10].
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 wt strain or the different mut strains. Consequently, the number of commensal species is more affected by the host immune attack than the number of Ye mut strains, which are, of course, more affected than the wt 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 system of differential equations [10]: 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, wt Yersinia, and mut 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 describes the growth rates of the bacterial populations, which are limited through a given capacity "C" and are defined as follows: C L mut growth rate in the lumen (13) Model parameters are shown in Table 2. Since the system models the Yersinia population in the mucosa and lumen, we assume that all state variables and parameters of the model are non-negative ∀t ≥ 0. The basic reproduction ratio R 0 is calculated by the fraction of the transmission rate for each strain (spill-over) and the average infectious period for each strain in compartments, as defined by Diekmann et al. [28].
In addition, van den Driessche and Watmough [29,30] 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. Suppose there are n disease compartments and m non-disease compartments, and let x ∈ R n and y ∈ R m be the subpopulations in each of these compartments. In this splitting, the authors thus represented the diseased compartment with F i (x) as the rate of the appearance of new infections in compartment i, and V i as the rate disease progression, death, and recovery reduced in compartment i, which can be written asẋ To understand compartment V i better, this can be written as V i = V − i + V − i , where V + i is defined as the rate of transfer of individuals into compartment i by all other means, and V − i is the rate of transfer of individuals out of the compartment i. Therefore, the disease transmission model consisting of non-negative initial conditions will be as follows: d dt It is assumed that functions F i (x) and V i = V − i + V − i are continuously differentiable at least twice in each variable, and they satisfy the assumptions A(1)-A(5) described below [29,31]: A(1) F i (0, y) = 0, V i (0, y) = 0 : ∀y > 0 and i = 1, . . . , n (no immigration of individuals into the disease compartments) A(2) F i (x, y) ≥ 0 : ∀x i ≥ 0 ∧ y i ≥ 0 and i = 1, . . . , n (the new infections will be represented by F , so it cannot be negative) A(3) V i (x, y) ≤ 0 : whenever x i = 0, i = 1, . . . , n (if the compartment is empty, it can only have inflow, and the net outflow from the compartment must be negative) The systemẏ = G(0, y) has a unique asymptotically stable equilibrium, y * (all solutions with initial conditions of the form (0, y) approach a point (0, y * ) as t → ∞) Proof. See [29].
denotes the spectral radius of a matrix A. 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 matrix of the system. If the eigenvalue is equal to zero, the bifurcation is a steady-state bifurcation. Therefore, we now recall the analysis of the center manifold near the critical (X = X 0 , R 0 = 1), which allows clarifying the direction of the bifurcation near the bifurcation point using the Center Manifold Theorem 2. This theory describes the local stability at the non-hyperbolic equilibrium (linearization matrix has at least one eigenvalue with zero real parts) and the existence of another equilibrium (bifurcated from the non-hyperbolic equilibrium).

Center Manifold
The center manifold theorem provides a means for systematically reducing the dimensions of the state spaces, which need to be considered when analyzing bifurcations of a given type.

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 , andσ u , with Let the (generalized) eigenspaces of σ s , σ c , and σ u be E s , E c , and E u , respectively. Then, there exists C r stable and unstable invariant manifolds W u and W s tangent to E u and E s at 0 and a C r−1 center manifold W c tangent to E c at 0. The manifolds W u , W s , and W c are all invariant for the flow of f . Stable and unstable manifolds are unique, but W c need not be.
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 manifold and then reduce the system to its normal form. Without loss of generality, we assume that x = 0 is the fixed point of interest for the system. Suppose J has n c eigenvalues with zero real-parts and n s eigenvalues with negative and positive real parts, 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 parts, the Center Manifold Theorem 2 guarantees that there exists a smooth manifold 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 (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( (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 if the type of transcritical bifurcation is forward or backward. Proposition 1. Assume that conditions A(1)-A(5) are satisfied. Furthermore, assume that the following hypotheses are satisfied by the system in (14) and (15) Then, the transcritical bifurcation of the system in (14) and (15) at R 0 = 1 is forward, and if at least one of these features is not present in the model structure, then a backward bifurcation may occur [31].

Results
We used the model in (1)-(7) to understand how the dynamics change following variations of different parameters. We calculated the equilibria of (1)-(7), conducted a linear stability analysis, and identified the analytical conditions that lead to a transcritical bifurcation.

Existence of Equilibria
For mathematical convenience, we divide the model in (1)-(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 Ye growth rates α (wt) and α (mut) , respectively. 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 (1)-(7) would be as follows: • The trivial equilibrium point is as an origin equilibrium (0, 0, 0, 0, 0, 0, 0). This solution appears when all populations are extinct. For all parameters, this point never becomes stable due to the positivity of eigenvalues in (A2). • The first equilibrium point appears in the absence of Yersinia Y (1)-(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 disease-free state P 0 to be biologically meaningful, the conditions γ < α (B) and β < α (B) must hold. These conditions correspond to the maximal growth rate of intestinal bacteria exceeding the rate at which intestines are charged and the maximal immunity action, which is not that strong in the absence of Yersinia strains. However, the population of the immune system is at its maximum carrying capacity (in health, not in fighting with any infection). • A second equilibrium corresponds to the commensal bacteria's persistence and the Yersinia mut strain in the absence of the wt strain. Without loss of generality, the commensal bacteria are supposed to be zero because they are not infective. This point is obtained by setting Y • The other equilibrium corresponds to the persistence of commensal bacteria and the Yersinia wt strain in the absence of the mut strain. Without loss of generality, the commensal bacteria are supposed to be zero because they are not infective. This point is obtained by setting Y • Finally, the last equilibrium point corresponds to a state of the co-existence of wt and mut Yersinia strains. This point is achieved by supposing B M = B L = 0: where Z is defined as 1 − to make the equilibrium point biologically meaningful.

Analysis of the Disease-Free Equilibrium Point
The system's behavior close to the equilibrium points was determined through linear stability analysis and bifurcations. It was carried out by calculating the Jacobian matrix of the model in equilibrium points. The Jacobian matrix for (1)-(7) is given by Equation (A1) in the Appendix A.
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 describing an infectious disease transmission within a heterogeneous population, the host population is grouped into two general classes: the infected and uninfected compartments. Therefore, the system (1)- (7) is divided into two infected and uninfected compartments.
In classical epidemic models, it is common to observe that a disease-free equilibrium loses its stability for R 0 = 1, and a transcritical bifurcation occurs. A transcritical bifurcation is when a fixed point exists for all parameter values and is never destroyed. However, as the parameter values vary, this fixed point transitions from a stability region to an instability region. Biologically speaking, as R 0 < 1, the disease-free equilibrium would stay stable. Thus, only the commensal bacteria persist, and wt and mut Yersinia strains cannot pass the invasion boundary. Therefore, as R 0 increases, wt and mut Yersinia strains, fed by commensal bacteria, can appear. We mathematically analyze this aspect within the structure of the model to assess which parts of the structure might be responsible for the occurrence of the transcritical bifurcation. Let us consider system (1)-(7) with the above-calculated eigenvalues. Using the eigenvectors (A9) shown in the Appendix A as a basis for a new coordinate system, we set the transformation matrix whose columns are the eigenvectors; Under this transformation, the system (1)-(7) changes to 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}, satisfying the following equations: Thus, the reduced system isu The advantages of a center manifold are clear from this calculation. We may study a one-dimensional system instead of a seven-dimensional system. That is, the method of center manifolds enables one to reduce the dimensions of the system by studying the flow restricted to the center manifold, in which the transients associated with the nonzero eigenvalues have decayed. As long as α (wt) < α (B) and α (mut) < α (B) , the disease-free equilibrium point will persist. When the infection rate of wt and mut strains increases, the system becomes unstable at the disease-free equilibrium point, and the trajectory of the system approaches asymptotically to the bacterial population of wt or mut and, finally, to the co-existence of the strains.

Computing and Analysis of the System through the Basic Reproduction Number
For scenarios where multiple strains (subtypes) of an infectious disease exist, an R 0 number is calculated for each strain. Hence, in our model, we must define two different R 0 , one for the wt and the other for the mut strain, rather than an R 0 for the whole model. This is a significantly more difficult task. Herein, we introduce the following variations to the basic reproduction number, which are calculated by the average number of secondary infections: To understand the role of the basic reproduction number, we also define a single reproduction number for commensal bacteria R com 0 , computed as an expected number of secondary cases of commensal bacteria produced by a single commensal bacteria: Due to biologically meaningful disease-free state P 0 and holding the conditions β < α (B) and γ < α (B) , the disease-free state would be meaningful if R com 0 > 1. Therefore, the commensal bacteria appear, R com 0 slightly reaches above one, and a stable disease-free state happens. Due to the presence of commensal bacteria and instability of the trivial solution, the solution corresponds to the extinction of all populations that never appear because if R com 0 < 1, all components of the disease-free equilibrium will lie out of the biologically meaningful region.
For R wt 0 < 1 and R mut 0 < 1, thus, λ 1 = 0 is a simple zero eigenvalue, and the other eigenvalues are real and negative. This implies that transcritical bifurcation occurs in the disease-free equilibrium 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 resulting in the Jacobian matrix having a single eigenvalue that is equal to zero, then either the mut strain or the wt appears and persists for R mut 0 or R wt 0 , respectively. As shown in Figure 2, when the disease-free equilibrium loses its stability, different scenarios can occur. Herein, we discuss the different cases. . Thus, the intersection of the transcritical curves R wt 0 and R mut 0 results in a triple transcritical bifurcation. As shown in (A9), the Jacobian has a triple zero eigenvalue at this point (λ 2 = 0, λ 3 = 0). Kuznetsov [34] has proved that such a point would be an indicator of the onset of a non-degenerate or degenerate Bogdanov-Takens bifurcation [34,35]. The disease-free equilibrium P 0 loses its stability, and the wt-free and mut-free include one simple zero eigenvalue (λ 3 = 0), meaning that the dynamics of the model change as the target parameter is within the threshold value. (II) If R wt 0 > 1, the wt strain equilibrium in region II will persist when R mut 0 < 1. The wt strain will spread and possibly persist within the host population. In general, for a strain to persist, its basic reproduction number has to be strictly greater than one. Therefore, in this region, the disease-free, mut strain, and co-existence state exchange stability: P 0 becomes unstable, P 2 becomes locally asymptotically stable, and P 1 and P 3 remain unstable. This means that the immune system could kill one of the strains more efficiently. (III) If R mut 0 > 1, the mut strain equilibrium in region III will persist when R wt 0 < 1. The mut strain will spread and possibly persist within the host population since its basic reproduction number is greater than one. Therefore, in this region, the disease-free, wt strain, and co-existence state exchange stability: P 0 becomes unstable, P 1 becomes locally asymptotically stable, and P 2 and P 3 remain unstable. This means that the immune system could defeat the wt strains. However, the risk of this situation to happen is low because the mut strains are influenced more efficiently than wt strains by immune action. (IV) If R wt 0 > 1 and R mut 0 > 1, the co-existence population spreads, and both strains persist. The overall R 0 can be defined as R c while one with d = 0 is completely recessive; scenarios of incomplete dominance (0 < d < 1), under-dominance (d < 0), and over-dominance (d > 1) are possible as well. For instance, a mut could achieve a higher R 0 than the wt via a higher growth rate that increases transmission. In a co-infection, the faster-growing mut strain would outcompete the wt and reach its maximum capacity. This situation would change the co-infection to the conditions where a single infection happens. Thus, the overall R 0 of the co-infection would be similar to that of the mut by itself, making the mut a dominant one. Furthermore, the effort for having a co-existence equilibrium and analysis of the co-infection model will fail. By contrast, let us assume a mut strain achieves a higher R 0 . Nevertheless, the virulence of the wt strain neutralizes the higher R 0 value of the mut. This would make the mut a recessive one. In summary, virtually any two-strain co-infection model can be mapped to a set of values for d, allowing scenarios of particular interest to be explored in a context broader than the one possible with typical models. The mathematical model (1)-(7) calculates the number of pathogen-specific characteristics in different layouts, e.g., when colonization resistance mediated by the endogenous microbiome is lacking or when the immune response is partially impaired. In this paper, we use the experimental data obtained in the lab upon infection of a host either lacking a microbiome (mimicking antibiotic treatment of patients) or a fully functional immune system [10].
To challenge the numerical simulation with experimental data, we tested whether the numerical simulation from the mathematical analysis could fit the experimental data.
Correspondingly, we simulated the process of Yersinia enterocolitica co-infection in specific-pathogen-free (SPF) (i.e., wt), germ-free (GF), and MyD88-deficient (MyD88 −/− ) mut mice. Using the model values in Table 3, which was estimated through experimental data [10], we first examined the influence of the infection rates α (B) , α (wt) , and α (mut) on the dynamics beginning by constructing one-dimensional bifurcation diagrams using α (B) as the bifurcation parameter and fixing α (wt) and α (mut) . The model values in Table 3 are the output of estimation parameters in Table 2 of the (1)-(7) by experimental data [10]. In theory, some parameter values were defined based on the biological aspect, and the rest were estimated by the optimization problem with the maximum likelihood estimation [10].
As shown in Figure 3, R com 0 should be larger than one to achieve the disease-free equilibrium. As long as R wt 0 and R mut 0 are smaller than one, the disease-free equilibrium is stable.    Table 3 in model (14), the commensal bacteria in the mucosa appear when R com 0 > 1. Therefore, as long as R com 0 < 1, only the trivial solution for the model exists. Since the trivial solution is always unstable, the extinction of all populations is never achieved. (b) By fixing the parameter values Table 3 in model (14), the commensal bacteria in the lumen appear when R com 0 > 1. Therefore, as long as R com 0 < 1, only the trivial solution for the model exists. Since the trivial solution is always unstable, the extinction of all populations is never achieved.
Following the analysis of the model with respect to α (B) , two cases are considered. The first case is when α (wt) and α (mut) are equal. In this case, we do not have co-existence of wt and mut (the denominators of Y 3(wt) L and Y 3(mut) L will be equal to zero). The second one is when α (wt) and α (mut) have different values as the co-existence equilibrium is achieved and is biologically meaningful.
Thus, we start with values of the wt and mut-strain infections rate given by α (wt) = 4.44 × 10 −1 and α (wt) = 4.44 × 10 −1 in Table 3. Corresponding to those infections rates, the basic reproduction numbers R 0 (wt) and R 0 (mut) in terms of α (B) are computed.
Through the definitions of R com 0 , R wt 0 , and R mut 0 , we display an x-scale as α (B) in Figure 4. This shows the appearance of a different equilibrium of Model (14) in terms of the basic reproduction numbers.
As shown in Figure 4, we expect no co-existence equilibrium when α (wt) = α (mut) , and only wt equilibrium exists (mut-free equilibrium) when 1.76 < α (B) < 11.20. However, the disease-free equilibrium can happen when α (B) is large enough to surpass α (wt) and α (mut) . These findings are displayed in Figure 5.
The appearance of a freedisease equilibrium point In this case, the region behind the blue point is not biologically meaningful. We do not have an equilibrium point that corresponds to a state of co-existence of wild-type and mutant. The condition is necessary for the existence of a co-existence equilibrium.
( α (wt) ≠ α (mut) ) The appearance of the equilibrium point corresponds to a state of co-existence of wild-type and mutant Yersinia strains.
This region is not biologically meaningful. The only equilibrium in this region is the trivial solution.
The appearance of equilibrium corresponds to the persistence of the Yersinia wild-type strain in the mutant strainsʼ absence.  Table 3 in Model (14).  The sensitivity of parameter α (B) with respect to commensal bacteria, wild-type, and mutant Yersinia in the mucosa and lumen for 336 h. 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 a wt equilibrium (R wt 0 > 1). Thus, (a,e) are increasing fast and stay at the maximum level as (a,d) are going back to zero. Additionally, (c,d) do not grow anymore (R mut 0 < 1).
Secondly, to show the appearance of the co-existence, we consider the parameter values of in Ye wt/T3S0 from Table 3 with the assumption of 1.86 = α wt = α mut = 1.83. Thus, we compute R wt 0 and R mut 0 with respect to parameter α (B) . This results in the appearance co-existence when 1 < α (B) < 4.90 as Figure 6 shows.  Furthermore, we analyze the immune system's influence. In our model, at least one Ye within the mucosal compartment activates the immune system. This activation increases proportionally to the number of Ye. The immune system is assumed to influence the bacterial populations primarily at the mucosal site compared to bacteria within the lumen. Herein, we simulate the process of Yersinia enterocolitica with an immune response that was experienced in SPF, GF, and MyD88 −/− mice. As the responses of the immune system in these mice are different, we conclude different behaviors. Let us denote the numerical simulation: • The effect of the maximum rate of immune growth κ on wt Yersinia strain in the mucosa, Figure 7; • The effect of the maximum rate of immune growth κ on mut Yersinia strain in the mucosa, Figure 8; • The effect of the maximum rate of immune growth κ on wt Yersinia strain in the lumen, Figure 9; • The effect of the maximum rate of immune growth κ on mut Yersinia strain in the lumen, Figure 10, The experimental data obtained in the lab showed that the Ye mut strain was eliminated more efficiently than Ye wt in the mucosal compartment. However, since GF mice lack a microbiome, both the Ye wt and the mut strains can colonize the gastrointestinal tract (GIT) at high numbers. In the SPF-colonized MyD88 −/− mice, the immune response is weaker [10]. Therefore, as shown in Figures 7-10, the following results conclude: • Figures 7d and 8d show a slow reduction in the wt strain in comparison to mut strain as κ changes. • Figures 7b, 8b, 9b and 10b in comparison with similar situations in SPF and MyD88 −/− show elimination of wt and mut strains is less efficient. • Figures 7c, 8c, 9c and 10c show that wt and the mut strains increased faster. Besides, the influence of κ on wt and the mut strains cannot project properly (Figures 9f and 10f) at the same speed of producing strains. Our results allow us to state that if mice possess a healthy microbiota and immune system, as long as the growth rate of commensal bacteria is larger than the growth rate of wt and mut Yersinia, then the infection will not spread, as Ye strains cannot enter the lumen compartment. Therefore, the disease-free equilibrium, P 0 exists and is stable when R wt 0 and R mut 0 are less than one while α (B) is large enough. To challenge this situation, we tested the data values Ye SPF wt/A0 from Table 3 with the assumed boundary α (B) > 11.20 to face a disease-free state, Figure 4. Our computer simulations ran for 366 h and are shown in Figure 11. As predicted for the disease-free state, the commensal fractions B M and B L approached their carrying capacities C M and C L , respectively.
(a) Appearance of the disease-free state (b) Projection of the disease-free state Figure 11. Diagram displaying the disease-free state by fixing the parameter values Ye SPF wt/A0 from Table 3 in Model (14). As long as α (B) is large enough, the disease-free state persists. However, by reducing α (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 (14), as shown in Figure 4. Therefore, we note that the disease-free equilibrium always exists and is (asymptotically) stable whenever the reproductive numbers for both strains of the disease are below unity. Simulations support these findings. The two single-strain equilibria, where one strain persists while the other strain dies out, are symmetrical and exist when at least one of the strains has a reproductive number above unity. However, in our model, (14) with parameter values from Table 3, we could not observe a state that corresponds to Ye mut strain equilibrium. As shown in Figure 4, there is no region where R mut 0 > 1 and R wt 0 < 1. Whenever R mut 0 > 1, R wt 0 is already above one; this facilitates the situation for a co-existence equilibrium. This is where both strains persist, existing when both reproductive numbers are above a certain threshold. However, we could not analytically solve the stability criteria for the equilibrium due to complexity, but simulations show that these stability criteria exist.

Discussion and Conclusions
In this paper, we proved that our computational model of the Yersinia enterocolitica [10] population in co-infections of mice with Ye wt and mut strains always showed stable results in terms of R 0 , such that if we can keep or reduce R 0 s to less than one, then wt and mut strains cannot spread, and the infection resolves. The medical control strategies for this infectious disease, like other infectious diseases, are frequently predicted through the basic reproduction number. To know what the difference is, we should be aware that the decomposition in dynamics and the designation of compartments as infected or uninfected are unique, and the basic reproduction number is achieved in the same way for all different compartmental models. This also raises the question of how the infection becomes persistent or is resolved based on the calculation of R 0 as a significant measure in control policy.
The answer to this question depends on the magnitude of the basic reproductive number, R 0 , since calculating the stability of an equilibrium is very complicated. Therefore, we computed R 0 for our computational model based on different strains. This resulted in three R 0 s: one R com 0 , R wt 0 , and R mut 0 .
As it is obvious in any compartmental model, the disease-free equilibrium is locally (asymptotically) stable when 0 < R 0 < 1 and unstable if R 0 > 1. In other words, when 0 < R 0 < 1, every infectious strain will cause less than one secondary infection; hence, the strain will eliminate. When R 0 > 1, every infectious strain will cause more than one secondary infection; hence, Yersinia infection strains will persist. However, all public health control measures can usually be based on methods that, if practical, would lower R 0 to below unity. On the other hand, the co-infection equilibrium is locally stable when R 0 > 1 and unstable when 0 < R 0 < 1.
This trivial result is essential, but we had to adjust our multi-strain model with different R 0 s. To more smoothly analyze R 0 s, we looked for a common parameter in all R 0 s. This resulted in the parameter α (B) . We thus analyzed at which level of sensibility the system's parameter α (B) changes the model's invasion behavior. This came to an end when a disease-free state occurs, when α (B) satisfies the conditions β < α (B) and γ < α (B) . In addition, if α (B) is large enough, then the disease-free state can be stable, and the infection with other strains does not appear or spread. A mut strain equilibrium does not occur for parameter values estimated through experimental data [10] since there is no region with only R mut 0 > 1. Although there is a region with two threshold parameters in the model, R wt 0 and R mut 0 stay above one for a co-infection setup. This region includes a small range for the parameter α (B) . In addition, it is important that the maximum growth rates associated with the wt and mut strains should be unequal. Otherwise, the co-existence scenario cannot take place.
Furthermore, we looked for conditions in the immune system to see different scenarios of a weaker and more robust immune system in fighting against infection. This resulted in the analysis of systems with respect to parameter κ. We found that if mice possess a healthy microbiota and immune system, as long as the growth rate of commensal bacteria is larger than the growth rate of wt and mut Yersinia, then the infection will not spread, as Ye strains cannot enter the lumen compartment.
Our results indicate that the within-host dynamics of immunity can, in principle, have significant consequences on population-level dynamics. However, immunity alone never creates a backward bifurcation of the disease-free steady state under biologically realistic hypotheses, and this would require some other complementary conditions.
Altogether, Yersinia is a great model system and can be used to predict the spread of infections for more clinically relevant bugs, i.e., enteropathogenic or enterohemorrhagic Escherichia coli. To control the spread of infection, efforts must be made to ensure that the co-infection equilibrium is unstable, i.e., 0 < R 0 < 1, and in our model R wt 0 < 1 and R mut 0 < 1. Furthermore, an additional factor to take into account is whether antibiotics are used from the beginning of the infection, in which case the spread of Yersinia may or may not be extinguished. We also believe that further research into the elaborate makeup of the strains and how the immune system differentiates between strains could be incorporated to more accurately represent control policy.

. The Eigenvalues mut Strain Equilibrium
The eigenvalues of equilibrium of mut strain are calculated as Appendix A. 6

. The Basic Reproduction Numbers
The basic reproduction number is calculated as follows: 1. Based on Equations (14) and (15) 3. F i (x) corresponds to the rate of the appearance of new infections, and V i corresponds to the rate of transfer assumed as and 4. The basic reproduction number is achieved as the spectral radius of the matrix FV −1 ; , and x j are as follows: x j = Y  [26] (Level 3 Version 2 [27]) from BioModels database [36] under model identifier MODEL2002070001.