Epithelial stratification shapes infection dynamics

Infections of stratified epithelia contribute to a large group of common diseases, such as dermatological conditions and sexually transmitted diseases. To investigate how epithelial structure affects infection dynamics, we develop a general ecology-inspired model for stratified epithelia. Our model allows us to simulate infections, explore new hypotheses and estimate parameters that are difficult to measure with tissue cell cultures. We focus on two contrasting pathogens: Chlamydia trachomatis and Human papillomaviruses (HPV). Using cervicovaginal parameter estimates, we find that key infection symptoms can be explained by differential interactions with the layers, while clearance and pathogen burden appear to be bottom-up processes. Cell protective responses to infections (e.g. mucus trapping) generally lowered pathogen load but there were specific effects based on infection strategies. Our modeling approach opens new perspectives for 3D tissue culture experimental systems of infections and, more generally, for developing and testing hypotheses related to infections of stratified epithelia.


A. Analysis and parameter exploration.
A.1. Uninfected epithelium. The uninfected epithelium model is a linear system of ODEs that can be studied analytically. We use a continuous time model because the stages overlap in time and have different birth/death rates. The system has only one equilibrium, which is Its Jacobian matrix is a diagonal matrix, and thus the eigenvalues are easily found (λ1 = −µ, λ2 = −ν, λ3 = ∆q ρp).
Since µ, ν, and ρp are biological parameters and must be positive, this equilibrium is stable iff ∆q ρp < 0. For this condition to hold (and the epithelium size to make biological sense), ∆q must be negative. The fact that the eigenvalues are never complex indicate that there are no limit cycles and that the solution grows (or decays) unbounded when the equilibrium is unstable.
The eigenvalues also allow us to calculate the relaxation time, which indicates which parameters have the most effect on how quickly the system reaches equilibrium (i.e. how quickly the epithelium reaches homeostasis). If the eigenvalues are all negative, the leading eigenvalue, λ l , determines the long-term behaviour of the system since it has the largest relaxation time, −1/λ l . Which eigenvalue is the leading eigenvalue changes depending on which parameter values are used. Therefore, we used a sensitivity analysis to study within our given parameter space which parameter has the most effect on time to equilibrium. Our results show that ∆q is the most important parameter ( Table 2 in main text).
Next we investigate what governs the thickness, or in other words the sum of the equilibria. Let h be a measure for how thick the epithelium is at homeostasis (h =Ũp +Ũ d +Ũs).
To determine which parameters affect thickness the most, we maximize h. To do so, we simplify the model by making two assumptions: i) we do not follow Us because it simply collects cells leaving the U d population and ii) we re-parameterize the model, such that, C1 = ρ b (1 − ∆p)N b , C2 = ρp ∆q, and C3 = ρp(1 − ∆q). Together this give the following simplified model, Using this set of equations, we get h = C 3 C 1 C 2 ν + C 1 C 2 . By computing the partial derivative of h with respect to the 4 parameters and setting them to 0 (i.e. solving for ∂h ∂C 2 = ∂h ∂C 3 = ∂h ∂C 1 = ∂h ∂ν = 0), and then checking the sign of the second order derivatives for concavity, we find that there is no finite set of parameter values that maximizes h. Therefore, to study how much parameters affect the maximum of h in a restricted parameter space, we use uncertainty and sensitivity analysis (as described in the methods).

A.2. Parameter inference.
In order to fix the carrying capacity of the basal layer, N b , we choose the maximum value of the DAPI counts of the basal cells in FOVs across all replicates. Experiments that reach homeostasis within the experimental time window will best estimate this parameter. Given that ∆p was estimated to be ≈ 0, this implies that a growing equation for the basal layer (equation 2 in the main text) is not necessary and future fits could just use the original model that assume a constant population size of U b .
In order to reduce the number of parameters estimated, we used the counts of cells dividing at each time point from the BrDU stains to determine a relationship between the replication rates as a means. This constrained the value of ρp to be smaller than the replication rate of basal cells, ρ b . The BrdU stains give the proportion of cells dividing in each layer (basal, non-keratinized suprabasal and keratinized suprabasal) at each measured time point and we assumed that this is represents the lit-up cells during the cell cycle stages of DNA systhesis (S), postsynthesis (G2) and mitosis (M). Following (1) the total time of the mean cell cycle, Tc, is the sum of the time in these three phases, Tr, and the time in postmitotic phase (G1), i.e. Tc = Tr + Tg1. The proportion of cells that are labeled and lit up (which we have from the BrdU data) are Li, where i = b or p for basal cells or parabasal cells respectively. We assume then that Li = Tr/Tc i and that Tr is constant (only the G1 phase can vary in duration) and therefore is the same for both basal and parabasal cells. Assuming the replication rates follow these relationships with the mean cell cycle time, ρ b = ln(2)/Tc and ρp = ln(2)/Tc, and substituting Tc b = Tr/L b and Tc p = Tr/Lp we obtain the following relationship ρp = ( Given that the data does not distinguish parabasals from differentiated cells in the non-keratinized suprabasal layer, we correct this equation with the fraction of parabasal cells,

Up
Up+U d , so that ρp is not applied to differentiated cells. This gives the equation we used for parameter estimation, Murall et al.

of 8
where Up and U d are the variables parabasal and differentiated cells in Model 1 in the main text. The mean proportions of dividing cells for each cell type was calculated from the BrdU stains, giving Lp L b = 0.02 0.15 . Thus the basal layer had more dividing cells at any given time than parabasal layer. This relationship had a strong effect on the inferred ∆q value and thus it should be estimated from the data for specific systems in future studies.
The inferred death rate, µ, was found to be nearly zero (Table 1, main text) which is consistent with the fact that in the experiment there is no mucus to wash away dead cells thus this parameter would expected to be zero experimentally (the dead cells pile up and would be counted).
Finally, using a loglikelihood ratio test, we found that for this dataset, a reduced model where ∆p = 0 outperformed our full model, however, other reduced models where ∆q = 0 did not, signifying then that ∆p is not necessary but ∆q is to explain our data.

A.3. Human papillomaviruses. For low-risk HPV infections, we
show that re-seeding (via some infection of the basal cells during at least the initial part of the infection) is necessary for creating a wart-like excess number of cells. We make the infection rate a decaying function that approaches zero, i.e. β(t) = β0e −bt , thus mimicking rapid microabrasion repair. In Figure S1b, we see that this quick removal of re-seeding decreases the number of infected cells (in this case, nearly 2 orders of magnitude lower than the baseline). As a consequence, the infection goes nearly undetected by the immune system and lasts for many years. Thus, if microabrasions repair quickly then higher burst sizes (i.e. HPV types with a colonization strategy) would be more adaptive.
In Figure S2a, b and c, we demonstrate how increasing either burst size or HPV-driven proliferation of basal/parabasal cells leads to faster growing infections and thus detection and clearance by the immune system. To model progression Figure  S2, we made ∆q, ∆p, m and µ increase logistically with respect to time using the theta-logistic equation (2). Let N be a dummy variable, which can be replaced by ∆q, ∆p, m or µ, has the form where K = 1 and r = 0.01. In order for the curves reach saturation more slowly θ = 1/2. These new time variant equations were rescaled so that they increased within the following ranges: ∆q [-0.012,0], ∆p [-0.012,1], m [0.1,1] and µ [0.67,4]. In addition ρp was increased to 3. Note that this causes the parabasal cells to bulk-up. In any case, a full investigation into modelling progression of HR-HPV is beyond the scope of this paper. In particular it would require changing several assumptions over time and including a more dynamic interaction with the immune system. Figure S3 we see the effect of when infection rates are different between the stages (i.e., β b < βp < βu), namely that there is asynchrony of the oscillations between the uninfected cells and the infected cells. Comparing the dynamics of U d and I d or Up and Ip helps visualize the timelag between the oscillations. From a dynamical perspective, this echoes ecological predator-prey results (3): if a 'predator' (here chlamydia) feeds upon a stage-structured 'prey' population (here epithelium) with differing attack rates between the stages (here β b < βp < βu), dampened oscillatory population dynamics can be observed. The time-lags between the epithelial layers and the asynchrony between the uninfectecd and infected epithelial layers stabilize the infection dynamics. If the infection rates are similar for all stages (β b = βp = βu), the dynamics are less stable, with large oscillations around zero that generate infections with only acute phases.

A.4. Chlamydia. In
A.5. Non-stratified HPV infection model. In order to compare what added benefit the stratification of the epithelium can give when it is explicitly considered in our model, we build a non-stratified HPV infection model. As in 'conical' viral dynamics models (4), we consider four main populations: target cells (which in this case are epithelium cells lumped into one population), T , infected cells, I, free virus particles, V , and a population of effector cells from the adaptive immune response, A. Assuming that the epithelium is at homeostasis when the infection begins, we assume T is at a carrying capacity and thus constant. We considered two other formulations where we had an explicit T equation where N b is a carrying capacity and λ = µ. The results we present below did not change significantly when using either of these two T equations (not shown), thus we present here the simplest version of the model. The three equations of the system are [6c] The epithelium cells are infected by free virions at a rate β V T . HPV drives infected cells to proliferate, ρa I. Infected cells die at a natural rate of µ I and are killed by the interaction of immune cells, κ A I. Free virions are released when infected cells die naturally, µ θ I and are cleared by the mucus at rate ζ V . Finally, immune cells proliferate as a response to the density of infected cells, σ I A. The model schematic is below ( Figure S4). We maintained the parameter estimates the same as the ones used in the main text (Table 3), except we varied ρa and θ to study their effect on the dynamics of the infection. The main assumption contrasting wart-associated and HR types is that HR ρa is higher (due to their stronger cell transforming properties).
We find that the non-stratified model of HPV infection is not able to reproduce what is known about HPV infections, particularly when comparing wart-associated vs HR infections. HR-HPV types have stronger cell proliferation properties than wart-associated types, however, they produce flat infections that, on average, last longer. The non-stratified model gives a contradictory result: the HPV-proliferation rate, ρa has a strong effect on infection dynamics such that the stronger it is the more infected cells accumulate and quickly ( Figure S5  A and B-D). Therefore, the model erroneously predicts that HR-HPVs accumulate more cells and have shorter infections ( Figure S5