Modelling cell turnover within Air-liquid interface cultures of human airway epithelium

To determine and evaluate how cell-type specific infectivity might influence SARS-CoV-2 infection within HAE cultures, we developed a mathematical model that describes the cellular turnover and tissue differentiation within the HAE-culture system (Fig 1A). The model distinguishes between the three main cell types observed within human airway epithelium, including basal (B), secretory (S) and ciliated (C) cells. Basal cells continuously proliferate at a rate α dependent on the total cell density in the system (N = B+S+C). They follow a linear differentiation pathway into secretory and subsequently ciliated cells at rate λ_s and λ_c, respectively [20], with each cell type experiencing a cell type-specific death rate δ_x (Fig 1B, Materials and Methods). By fitting our model to experimental data on cell differentiation of ALI cultures of bronchial epithelial tissue [21], we estimated the rates of cellular turnover and differentiation, assuming that the system reaches confluence at around ~45 days after seeding (Fig 1C and Table 1). Given the complexity of the model and the available data, multiple parameter combinations are able to explain the observed dynamics with dependencies preventing individual parameter identifiability (S1 Fig). Nevertheless, we estimate consistently that basal cells have a longer half-life than ciliated and secretory cells, with the latter one having a ~3 (ciliated) to ~6 (basal)-fold shorter lifetime compared to the other two cell types (t_1/2b = 49.9 days vs. t_1/2s = 8.4 days vs. t_1/2c = 24.2 days). The observation that secretory cells have the shortest half-life of the considered cell types is in line with observations in vivo [20]. With our best parameterization, the bronchial epithelial tissue at confluence is dominated by ciliated cells, as generally observed [21], comprising ~46.9% [19.9%, 60.6%] of the cells, followed by basal (39.7% [28.4%, 60.2%]) and secretory cells (13.4% [9.4%, 22.0%]) (Numbers in squared brackets denote the 90%-prediction intervals).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Cell differentiation dynamics in ALI-culture systems of human airway epithelium.

(A) Sketch of human airway epithelial culture system. (B) Schematic of a mathematical model describing cell differentiation and turnover for human airway epithelial cell cultures distinguishing between basal (B), secretory (S) and ciliated (C) cells. For details see Materials and Methods. (C) Relative proportion of each cell type during differentiation of ALI-cultures of bronchial epithelium. Mean and standard deviation of the experimental data from [21] (dots/whiskers), with the predicted dynamics obtained when fitting the model shown in (B) to these data. The best fit (solid line) and 90%-prediction intervals for each cell type are shown (basal-green, secretory-blue, ciliated-orange). (D) Cellular Potts model of ALI-culture system mimicking pseudo-stratified epithelium at different time points during differentiation. Basal cells (green) are “overgrown” by secretory (blue) and ciliated (orange) cells. For graphical clarity, a reduced cell culture system is simulated starting with only ~1265 cells in total with N_max = 4085 cells. The average and range of n = 10 individual simulations using the best parameterization (Table 1) are shown with numbers indicating the average frequency of each cell type around day 80 post seeding (basal: 40.3% [39.3%, 42.5%], secretory: 13.6% [12.7%, 14.8%], ciliated: 46.1% [44.3%, 47.8%] (mean [min-max] across 10 simulations) (see also Movie1 at https://github.com/GrawLab/SARS-ALIculture).

https://doi.org/10.1371/journal.pcbi.1011356.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameter estimated for development and differentiation of bronchial epithelial ALI-culture systems.

Parameters were estimated by fitting the model for epithelial differentiation dynamics (Eqs (1) and (2), Fig 1B) to the data for untreated ALI-cultures [21]. Values indicate best estimate and 90%-credible intervals (CrI) of estimates using an Approximate Bayesian Computation approach (pyABC, [23]) (see also S1 Fig).

https://doi.org/10.1371/journal.pcbi.1011356.t001

Using the parameterization obtained above, we developed an individual cell-based model to account for the spatio-temporal dynamics of epithelial regeneration and cell differentiation. The model is based on a cellular Potts modeling (CPM) framework [22] which enables us to simulate inter- and intracellular dynamics on a local level by accounting for the polarized orientation of cells, i.e., basal and apical side of the epithelium, and varying morphologies of individual cells with increasing cell densities (Fig 1D). Our model is able to mimic the development of pseudo-stratified epithelium with basal cells slowly being overgrown by secretory and ciliated cells that will dominate the apical side of the culture (Fig 1D, see also Movie1 at https://github.com/GrawLab/SARS-ALIculture). Thus, we developed a model system that allows us to simulate and analyze the spatio-temporal dynamics of HAE cultures on a single-cell level.

Determining cell-type specific SARS-CoV-2 infection dynamics within bronchial human airway epithelium

To extend our framework for determining cell-type specific infection kinetics of SARS-CoV-2 among human airway epithelial cells and analyze their influence on the progression of infection, we followed a two-step approach. In a first step, we extended our cell population-based mathematical model (Fig 1B, Eq (1)) to additionally account for the infection dynamics (Fig 2A). Analogous to standard models of viral dynamics [24,25], we distinguish between uninfected (X), infected (X_I) and infectious cells (X_J). Each of the different cell types considered (basal, B, secretory, S, ciliated, C), is assumed to experience cell-type specific rates of cell infection β_x, of turning infectious and productively infected κ_x, viral production ρ_x, and cell death of productively infected cells δ_Ix, with x∈{B,S,C}. Infection occurs proportionally to the total viral load, V, which is cleared from the system at rate c. The model and the mathematical equations are explained in detail in Materials and Methods. As an extension of this standard model, and in accordance with previous observations [12,26], we considered a proportion of each of the considered cell types to become refractory to infection (Fig 2A). This accounts for the spatial separation of cells within the culture and within different cell layers (basal vs. secretory/ciliated) in the light of the well-mixed assumption of the model, as well as for cells that are protected from infection, e.g. by innate immune mechanisms [26] or low cilia expression [12].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Cell-type specific infection kinetics.

(A) Schematic of a viral dynamics model for SARS-CoV-2 infection within human airway epithelium considering uninfected (X), infected (X_I) and infectious/productively infected cells (X_J) for each of the different cell types X∈{B,S,C}. In addition, a population of cells being refractory to infection (X_R) is considered for each cell type to account for spatial separation and possible innate immune protection of cells (grey arrows). (B-D) Measured (dots) and predicted (lines/shaded area) for the viral load (B), relative proportion of infected cell types (C) and total cell count for each cell type (D) by fitting the model shown in (A) to data on SARS-CoV-2 infection of ALI-cultures of bronchial epithelium analyzed by scRNA-seq [27] are shown. Model predictions for the best estimate (solid line) and prediction bands (shaded area) based on the 90%-credible intervals of parameter estimates in the last ABC-generation (Table 3) are shown. (E-G) Predicted long-term oscillatory dynamics of the individual components. (H) Snapshots of the individual-cell based model simulating SARS-CoV-2 infection dynamics within bronchial epithelial tissue considering ~4.7×10⁴ cells in total (see also Movie2 at https://github.com/GrawLab/SARS-ALIculture). (I-J) Relative proportion of infected cells (I) and viral load dynamics (J) with the mean (solid line) and maximal range (shaded area) over 10 independent simulated cultures.

https://doi.org/10.1371/journal.pcbi.1011356.g002

To provide a first assessment of cell-type specific infection kinetics, the model was fitted to data published by Ravindra et al. [27] showing time-resolved measurements on viral load and single-cell sequencing analyses for ALI-culture systems of bronchial epithelial tissue infected with SARS-CoV-2 (Fig 2B–2D). In accordance with other studies [12,15], Ravindra et al. [27] found that ciliated cells represent the majority of infected cells, with ~26% of these cells being infected at 3 days post infection. While our standard model was not able to explain the experimental data as it would predict faster growth kinetics in the frequency of infected cells than observed experimentally when aiming to match viral load dynamics (not shown), the extended model with consideration of refractory cell populations provided an appropriate representation of the observed dynamics of the viral load (Fig 2B) and the relative proportions of cells being infected (Fig 2C). In accordance with the observation of ciliated cells being the main target of SARS-CoV2 infection [1,12,27], our model predicts a substantial, rapid loss of this cell type with only ~50% of them remaining at ~60 hours post infection (Fig 2D). Given tissue regeneration, the model would predict long-term oscillatory infection dynamics as they have been observed elsewhere [12,15] that would eventually saturate (Fig 2E–2G). Parameter estimates reveal an average duration of the viral eclipse phase of ~5h (2.9h, 20.0h) and an average lifetime of productively infected cells of ~12.9 h (6h, 28h) with no indicative differences between the individual cell types (Table 2). The largest differences between cell types are observed in the fraction of refractory cells, f, which is due to the model formulation. With the fraction of refractory cells, f, not being independently interpretable from the viral transmission rate, β, but the latter ones largely comparable between cell types (Table 2), infectivity is substantially reduced for basal cells compared to ciliated cells. While this could be due to the different spatial location of these cell types within the pseudo-stratified epithelium, the high fraction of refractory secretory cells compared to ciliated cells points towards a cell intrinsic difference in viral infectivity, with ciliated cells being more prone to infection.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Parameter estimates for cell-type specific infectivity.

Measurements on viral load and cell-type specific infection based on single-cell RNA sequencing data of SARS-CoV-2 infected ALI-cultures of human bronchial airway epithelium [27] were analyzed using the model shown in Eqs (3) and (4) and Fig 2A. The best parameter combination and the 90%-credible intervals (CrI) in the last generation of an Approximation Bayesian Computation approach (pyABC) after 50 generations are shown (*values from literature). Epithelial turnover dynamics were parameterized according to the best estimate shown in Table 1 (see also S2 Fig).

https://doi.org/10.1371/journal.pcbi.1011356.t002

To disentangle the contribution of spatial and cell-intrinsic effects on the infection dynamics, we incorporated the viral processes into our spatially-resolved model of the HAE culture system. Starting from the parameterization obtained by the ODE-approach (Table 2), we additionally adapted the transmission dynamics, as here no specific refractory populations are considered with the CPM explicitly accounting for the spatial separation of cells (see Materials and Methods for details on the adaptation process). Furthermore, in contrast to the ODE-approach which solely attributed infection spread to cell-free virus, the CPM particularly accounted for the spread of infection by direct cell-to-cell transmission within the model simulations, which has been found to be responsible for roughly 90% of the infections within tissues [28] (see also S2 Table). We used a step-wise approach considering different model assumptions to test their appropriateness in simultaneously explaining different observations on SARS-CoV-2 infection dynamics within tissues that comprise longitudinal measurements on specific cell-type frequencies and viral load dynamics [12,27] (see Materials & Methods for a detailed explanation). As a first result, we found that the spatial separation of cells was not sufficient to explain the different infection dynamics of ciliated and secretory, as well as basal cells. In addition, the spatially explicit model indicated the necessity of a substantially higher infectivity of ciliated in comparison to secretory and basal cells to match experimental observations. Furthermore, to explain the early rapid increase of the proportion of infected cells [27] in combination with the not fully destroyed tissue integrity at later time points [12], inclusion of a rate limiting process was needed to improve explaining the dynamics, as e.g. representative for an immune response [26]. Our best model parameterization (Fig 2H), for which we focused on the cell-type specific infection dynamics (day 0–7 post infection), led to a reasonable representation of the dynamics for the proportion of infected cells per cell type (Fig 2I), as well as viral loads within the observed ranges (Fig 2J), by slightly underestimating the dynamics within basal cells, as well as the peak viral load. Based on previous estimates indicating that spread of SARS-CoV-2 relies predominantly on cell-to-cell transmission (~90% of infections [28]), we obtain local infection spread that, in combination with the regenerative capacity of the tissue, determine the oscillatory dynamics (Fig 2H, Movie2 at https://github.com/GrawLab/SARS-ALIculture and compare Fig S7 in [12]). Assuming infection spread to be more dominated by cell-free virus transmission, our model would predict higher peak viral loads and more substantial loss of ciliated cells at early time points (S3 Fig). Thus, combining data from different sources, our analyses provided a first physiologically relevant spatially-resolved representation of SARS-CoV-2 infection within bronchial epithelial ALI-culture systems accounting for different cell type dynamics.

Spatial effects and tissue pathology

One of the hallmarks defining severity of SARS-CoV-2 infection within patients is the inferred pathology mediated to airway epithelial tissue in the lung [31]. Quantifying tissue pathology and relating it to functional loss, as well as evaluating the regenerative potential, are essential to predict disease progression. Within experimental tissue culture systems, tissue integrity is usually quantified by the transepithelial electrical resistance (TEER), determined by the electrical resistance between electrodes positioned at the apical and basal side of the culture [32]. However, these types of measurement cannot be obtained from patients impairing the ability to relate tissue integrity measures to pathology and possible disease progression. To overcome this issue, we developed a novel method that allows us to quantify tissue integrity based on image data (Fig 3A). The method analyses cell type specific membrane connections and relates them to TEER measurements, allowing us to quantify the resistance mediated by the present tissue pattern (see Materials and Methods for a detailed description of the method). Using TEER measurements of differentiating bronchial HAE culture systems [21] to parameterize resistance of cell type-specific tight junctions, our simulated culture systems show a rapid loss of tissue integrity after SARS-CoV-2 infection with only 70% (range 66–75%) of TEER remaining at day 4 post infection compared to day 0, comparable to experimental measurements (Fig 3B) [12]. The integrity is slowly repaired afterwards reaching almost homeostatic conditions within 25 days after infection, in correspondence to the reconstitution of the individual cell populations during tissue regeneration (Fig 3C). However, the inferred regeneration dynamics and oscillations seem to be slower than observed experimentally [12], indicating the possibility for altered cell proliferation and differentiation dynamics during acute tissue damage (Fig 1 and Table 1). In general, we do not observe a direct proportionality between the total number of cells and the tissue integrity measured by TEER, as the particular cell types and their connectivity influence these values.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Determining tissue integrity.

(A) Schematic of relating experimental measurements of transepithelial electrical resistance to image data based on individual cell tight junctions and spatial patterns. (B) TEER-dynamics for simulations of bronchial tissue (grey dots/red line) in comparison to previous experimental measurements from the literature (blue line) [12]. Data points indicate the mean (dots/red solid line) and min-max range (whiskers) across 10 independent simulation runs. (C) Corresponding dynamics of total and cell type specific cell numbers across the 10 simulations (mean–solid line, range–shaded area). Variations are minor with regard to the population sizes making the ranges hardly visible.

https://doi.org/10.1371/journal.pcbi.1011356.g003

Influence of cellular composition on tissue pathology and regeneration dynamics

The replication efficacy of SARS-CoV-2 is observed to vary dependent on the specific respiratory tissue, with nasal epithelium showing the highest permissiveness [33–35]. To analyze the effect of tissue composition and regeneration dynamics on the progression of SARS-CoV-2 infection, we additionally applied our model on epithelial differentiation (Fig 1A, Eq (1)) to HAE culture systems of nasal epithelium [36], and bronchial epithelial cells exposed to different concentrations of cigarette smoke extracts (CSE) during development, to determine cellular turnover under these conditions. Hereby, the latter condition simulated impaired epithelial restructuration due to smoke exposure [21,37,38].

Compared to bronchial epithelium, nasal epithelium is characterized by a larger fraction of ciliated cells, with basal cells only making up ~18% at steady state compared to ~40% in bronchial tissue (Fig 4A compare to Fig 1 and S3 Table) [36]. Our analysis indicates that these different steady states are due to faster differentiation rates of basal into secretory, and secretory into ciliated cells, which are estimated to be roughly 2-fold higher in nasal compared to bronchial tissue with all other parameters being comparable (Tables 1 and 3 and Fig 4B). In contrast to nasal and healthy bronchial epithelium, airway epithelium of smokers is usually characterized by basal cell hyperplasia [39–41], with ciliated cells being less frequent (Fig 4A and 4B). This is also seen in our estimates on the differentiation dynamics of bronchial epithelium exposed to cigarette smoke extracts by slightly increased loss rates for secretory and ciliated cells (Table 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Tissue composition and infection kinetics.

(A) Composition of bronchial and nasal airway epithelium, including bronchial epithelium exposed to cigarette smoke extract (CSE) at different concentrations [21,36]. (B) Relative proportion of each cell type during differentiation of ALI-cultures of the particular epithelium based on previous experimental data (mean/standard deviation = dots/whiskers) [21,36] with model predictions fitting the regeneration model shown in Fig 1B to the data. The best fit (solid line) and 90%-prediction intervals for each cell type are shown (basal-green, secretory-blue, ciliated-orange) with parameter estimates given in Table 3. (C) Simulated dynamics for the relative proportion of specific cell types infected, viral load, and total cell numbers during SARS-CoV-2 infection given different ALI-culture conditions. The mean and range (min-max) over 10 independent simulations for each condition are shown. (D) Corresponding snapshots of infection patterns at specific time points (see also Movie 3–5 at https://github.com/GrawLab/SARS-ALIculture). (E) Relative size of viral load and total cell population of the different tissue conditions using bronchial epithelium as a baseline. The ratio of the means across 10 simulations each are calculated. (F) TEER measurements across 5 (nasal, CSE) and 10 (bronchial) independent simulations for each condition relative to each TEER measurement at day 0.

https://doi.org/10.1371/journal.pcbi.1011356.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Parameter estimated for development and differentiation of HAE systems for different conditions: Parameters were estimated by fitting the model for epithelial differentiation dynamics (Eqs (1)–(2), Fig 1B) to data from ALI-cultures of nasal epithelium [36] and of bronchial epithelial tissue exposed to 2.5% and 5% cigarette smoke extract (CSE)[21].

Values indicate best estimate and 90%-credible intervals (CrI) of estimates using an Approximate Bayesian Computation approach (pyABC, [23]) running for 15 generations with 2000 particles each.

https://doi.org/10.1371/journal.pcbi.1011356.t003

To determine how the different tissue compositions affect infection dynamics, we simulated SARS-CoV-2 infection within these conditions using the same infection parameters as obtained for bronchial epithelium (Tables 2 and S1), in combination with the tissue-specific regeneration and differentiation dynamics (Table 3). As expected, given the high availability of ciliated cells, infection within nasal epithelium is faster and more severe than in bronchial epithelium, leading to ~1.5-fold higher fraction of infected ciliated cells compared to bronchial epithelium at day 5 post infection, and a loss of ~80% of ciliated cells at 10 days post infection (Fig 4C, compare Fig 2I). With tissue regeneration being ~2-fold faster than in bronchial epithelium, higher viral load levels are maintained over longer periods of time supporting oscillatory dynamics and regrowth of infection, including continuous tissue reconstruction (Fig 4D). In contrast, bronchial epithelium exposed to cigarette smoke extract does not lead to more severe infection patterns with viral loads decreasing faster than in unexposed tissue (Fig 4C–4E, compare Fig 2I, 2J). However, tissue regeneration is partly impaired under these conditions, showing reduced ability of the tissue to rebuild homeostatic cell numbers after infection compared to untreated bronchial tissue (Fig 4E). This is not visible within the progression of the calculated relative transepithelial electrical resistance values, which only show a slightly slower progression rate of TEER normalization for CSE-treated compared to untreated bronchial tissue, while the faster regeneration capacity of nasal epithelium is confirmed (Fig 4F). In summary, our analysis indicates how tissue and cell type specific effects could affect SARS-CoV-2 infection dynamics within human airway epithelial tissue.

Air-liquid interface cultures

Nasal and bronchial ALI cultures were purchased from Epithelix (https://www.epithelix.com) and were received when they had been differentiated for 6 weeks and displayed the correct TEER. Once received the samples were placed in fresh media and allowed to recover. Two days post-delivery samples were washed to remove mucus from the apical membrane. The number of cells in an ALI culture was determined by fixing the transwells and staining for DAPI. The cell number was then generated by imaging on a Zeiss CD7 Cell Discoverer and nuclei counting with Illastik (see S1 Data and S4 Fig).

Model of cellular turnover and epithelial regeneration

To model the differentiation and cellular turnover within the Air-Liquid Interface culture system we consider the three main cell types of human bronchial epithelium including basal (B), secretory (S) and ciliated (C) cells. Basal cells continuously proliferate at a density-dependent proliferation rate, with a maximal proliferation rate α. The cells follow a linear differentiation pathway into secretory and ciliated cells at rate λ_s and λ_c, respectively [20]. All cell types die at a corresponding death rate δ_i with i∈{B,S,C}. The whole dynamics is then described by the following system of ordinary differential equations.

(1)

Hereby, N = B+S+C defines the total cell concentration and N_max the maximal carrying capacity of the culture system. To account for a Γ-distributed differentiation process, the system in Eq (1) was then extended, assuming that basal and secretory cells pass through k+1 intermediate steps before becoming a secretory or ciliated cell, respectively. The standard model for cell regeneration and differentiation is then defined by (2) with the total number of cells, N, calculated accordingly, i.e., N = B₀+B₁+⋯. For the analysis and parameter estimation, the number of intermediate steps was set to k = 4, i.e., having 5 different compartments (0,1,..,4), because it fits the data reasonably well and is a tradeoff between a non-exponential distribution (k = 1) and a computational low cost model to allow for reasonable distributions of differentiation durations.

Model of cell-type specific infectivity and viral dynamics

To model the infection dynamics of SARS-CoV-2 within the human epithelial cultures, the standard model for cell regeneration and differentiation (Eq (2)), is extended by the infection dynamics for each of the different cell populations. Basal, secretory and ciliated cells get infected proportional to the virus concentration V at cell-type specific infection rates β. After a viral eclipse phase with an average duration of 1/κ, infected cells I will turn into productively infected cells J and start exporting virus at rate ρ, having an average productive lifetime of 1/δ_I. For simplicity and given the fast turnover of infected cells, we neglect the loss of non-productively infected cells due to cell death. Each of the different rates of transmission, cell death, eclipse phase and viral production are initially considered to be cell-type specific. The additional equations to Eqs (1)/(2) (without accounting for the compartmentalization in the notation) are then given by (3)

Hereby, the parameter c defines the viral clearance rate. In addition, N is now defined by , i.e., determining the total concentration of cells within the system.

Analogously to cell differentiation, we also used compartmentalization for the viral eclipse phase, 1/κ, and the average lifetime of productively infected cells, 1/δ_I, with k = 4 intermediate steps each used for parameter fitting.

With this standard model for infection (Eq (3)) not able to explain the data (see main text), we extended the model by additionally considering a proportion of cells that is refractory to infection, X_R, X∈{B,S,C} for each of the three cell types considered. These populations account for the spatial separation of cells within the culture and within different layers (basal vs. secretory/ciliated), and could also be interpreted as cell protected by innate immune mechanisms [26]. The corresponding equations in Eq (3) above are then modified and extended by (4) with f_X, X∈{B,S,C}, denoting the fractions of the corresponding cell types that become refractory, and throughout. Please note that the compartmentalization for the processes of cell differentiation, the viral eclipse phase, and the lifetime of productively infected cells (analogously to Eq (2)) were applied as well, with Eq (4) then representing the model that best explained the experimental data.

Parameter estimation for cellular turnover and epithelial regeneration

We fitted the model shown in Eq (2) to the data of Schamberger et al. [21] (and de Borja Callejas et al. [36]), using an Approximate Bayesian Computation-Sequential Monte Carlo (ABC-SMC) method [23]. The model was fitted to the data that were obtained by digitalizing the corresponding figures considering the experimentally observed errors with an additional 2% Gaussian noise to correct for data points where errors could not be observed. In addition to the cellular frequencies, we also included a penalty term that ensured that absolute cell numbers maintained reasonable values when fitting the relative frequencies. As ALI-cultures reach confluence and full differentiated epithelium in between 28–50 days after seeding (e.g. https://www.epithelix.com/products/mucilair and https://www.stemcell.com/), additional penalty terms were added to ensure that the individual cell populations reached steady state at around day ~45–50. The cost function for each cell type is then defined by with t_i = 7, 14, 21, and 28 days defining the time points of measurements with actual observations, and μ_sim(t), μ_exp(t) and σ_exp(t) the predicted and measured mean and standard errors of the particular cell type, respectively. The parameters w₁ and w₂ define weights to enforce the system to reach a steady state at ~45 days. Here, we set w₁ = 0.025 and w₂ = 0.005 to ensure a smooth transition from the dynamics inferred from the experimental measurements (first term) towards the assumed steady state of the system based on qualitative assessments in the absence of comparable experimental data. The choice of the weights will impact the transition dynamics and, thus, the estimated model parameters: Larger values of w₁and w₂ will delay the model to reach steady state.

The complete cost function is then given by

The second term comprises the penalty for the total cell count, with N(45) giving the total relative cell count at day 45, and w₃ the weight set to w₃ = 0.2.

The cost function was minimized by Approximate Bayesian Computation (pyABC, [23]), using 5000 particles for initial calibration and 2000 particles per generation. Each evaluation run for up to 15 generations until sufficient convergence was reached. As expected, given the complexity of the model and the available data, correlation analysis indicates that not all parameters are identifiable, with especially the cell proliferation rate α and the maximal cell number N_max experiencing strong correlations (S1 Fig). Prior distributions for individual parameters were defined by α~U(0.1, 0.6), λ_s~U(0.05, 0.3), λ_c~U(0.05, 0.3), δ_b~U(0.001, 0.021), δ_s~U(0.02, 0.16), and δ_c~U(0.02, 0.07) with all units given per day.

Parameter estimation to determine cell-type specific infectivity

Parameters describing cell-type specific infectivity were estimated by fitting the model shown in Eqs (3/4) to the data obtained from Ravindra et al. [27]. These data contained information from scRNAseq on the frequency of infected basal, secretory and ciliated cells at 1, 2 and 3 days post infection of a bronchial HAE system with SARS-CoV-2.

The model was fitted to the data by Approximate Bayesian Computation (pyABC, [23]), running for 30 generations with 1000 calibration particles and 1000 particles per generation. The previously determined best estimates for bronchial epithelium were considered for epithelial cell turnover (Table 1). According to the measurements from Ravindra et al. at 1h post infection [27], an initial infection dose of V₀ = 2.8×10⁶ RNA copies was used.

A cellular Potts model of the ALI-culture system

We used the cellular Potts modelling framework within the software tool Morpheus [22] to simulate the spatio-temporal dynamics of infection and tissue regeneration within the pseudo-stratified epithelium of the HAE cultures system. The cellular Potts model is a grid-based, individual cell-based model which allows to simultaneously account for inter- and intracellular dynamics across multiple scales. Each cell consists of several grid sites and aims at reaching a specified target volume and perimeter. Heterogeneous cell morphologies are reached as a cell aims at minimizing its global energy, which is defined by its surface tensions and characterized by the so-called Hamiltonian [16,22]. The CPM implementation with its flexible cell shape allows us to account for varying cell volumes with increasing cell densities. This is determined by the parameters volume strength and the so-called aspherity, V_S and Ψ, which define the strength with which a cell tries to reach/maintain the desired target volume and surface area, respectively, and penalizes any deviations from them (see Morpheus and the documentation therein for a detailed explanation of the implementation [22]).

The model distinguishes between basal, secretory and ciliated cells. To simulate a pseudostratified epithelium with basal cells located at the basal membrane of the culture, V_S for basal cells is only 1/10 of the corresponding value for secretory and ciliated cells (S1 Table). With this parameterization, basal cells are “compressed” as secretory and ciliated cells increase in numbers, simulating the “overgrowth” of basal cells by the other cell types (see Movie1 at https://github.com/GrawLab/SARS-ALIculture).

Cell proliferation, differentiation and turnover

The dynamics of the individual cell types with regard to cell proliferation, differentiation and cellular turnover follows the individual processes already defined within Eqs (1/2) and Eqs (3/4) (see also Figs 1B and 2A). Parameterizations are based on the previously estimated parameters (Table 1) with stochastic simulations accounting for cellular heterogeneity. To this end, the different continuous rates were translated into probabilities for the discrete events within the CPM dependent on the used timestep size assuming a constant hazard model. For example, with δ_S defining the loss rate of secretory cells within the ODE-model, the probability of a secretory cell to die (i.e., to change compartments according to Eq (2) to account for a Γ-distributed lifetime) within a time step Δt in the CPM is calculated by . Upon cell death, dead cells are removed immediately, i.e., not inhibiting cell proliferation dynamics or the calculation of the transepithelial electrical resistance (TEER) (see below). Due to the dynamic cell morphologies within the CPM and to account for the pseudostratified tissue when modelling density dependent cell growth, proliferation of basal cells is influenced by the total number of cells within the cell culture, and not solely the local surrounding showing reasonable tissue dynamics and growth kinetics.

Estimating transepithelial electrical resistance (TEER) based on image data

To measure tissue pathology within the ALI-culture system we quantified transepithelial electrical resistance (TEER) in our simulations, which is a standard method to measure the integrity of a tissue [32]. It determines the electrical resistance of a tissue layer placed on a semi-permeable membrane by measuring the resistance of the cell layer between electrodes placed at the apical and basolateral compartments of the ALI-culture system. Assuming that the resistance per area of the cellular membranes is much larger than the resistance at the tight junctions (TJ) of the cells, R_membrane≫R_TJ, the total resistance of the tissue, R_Tissue, is the harmonic mean of the resistance of the tight junctions, since they are connected in a parallel circuit. If all tight junctions have equal resistance, this would give us: with M_TJ denoting the area covered by the tight junctions. The reported TEER-value is then defined by TEER = R_Tissue×M_total in Ω×cm², with M_total being the total area of the tissue, i.e., M_total = M_TJ+M_cells (see also [32]).

With three cell types (basal, secretory, ciliated) within the culture system, as well as uncovered area (medium), the total resistance of the tissue is then given by (5) accounting for all possible tight junctions for each of the different cell types with each other and the corresponding area covered, with i denoting the number of all pixels. Hereby, R_bas,bas, R_cil,cil, etc. define the resistance of the particular tight junctions, and N_bas,bas, N_bas,sec, etc. the number of pixels with the corresponding connections. To determine the area of each tight-junction combination, we used the images generated by our CPM simulations. For each pixel characterizing a tight-junction (i.e., shown in black), we determined the type of pixels in the direct neighborhood (3×3 area) to classify the type of tight junction dependent on the type of cells in the surrounding. For example, if all pixels belong to either ciliated cells or represent other tight junctions, this would count towards a cil-cil area. Using a resolution of 2000×2000 pixels for our images (1 image pixel = 0.6×0.6 μm²), the borders of the tight junctions are usually 1–2 pixels in width.

To parameterize transepithelial resistance values in our simulations, we used TEER measurements for the differentiation of bronchial epithelium under homeostatic conditions ([21], Fig 1C). For simplification, we assumed that the cell types within the different layers of the pseudo-stratified epithelium (basal vs. secretory/ciliated) behave similar:

This assumption leaves R_bas,bas, R_bas,cil, R_sec,cil and R_medium as unknown parameters that need to be determined. In addition, a factor accounting for the influence of infectivity of one of the cells to lower the resistance of the particular tight junction was determined as well. Parameters were obtained by adapting the calculation of the simulated culture systems based on Eq (5) to the TEER measurements of untreated bronchial epithelium under homeostatic conditions [21].