key: cord-262524-ununcin0 authors: Bankhead, Armand; Mancini, Emiliano; Sims, Amy C.; Baric, Ralph S.; McWeeney, Shannon; Sloot, Peter M.A. title: A Simulation Framework to Investigate in vitro Viral Infection Dynamics date: 2011-12-31 journal: Procedia Computer Science DOI: 10.1016/j.procs.2011.04.195 sha: doc_id: 262524 cord_uid: ununcin0 Abstract Virus infection is a complex biological phenomenon for which in vitro experiments provide a uniquely concise view where data is often obtained from a single population of cells, under controlled environmental conditions. Nonetheless, data interpretation and real understanding of viral dynamics is still hampered by the sheer complexity of the various intertwined spatio-temporal processes. In this paper we present a tool to address these issues: a cellular automata model describing critical aspects of in vitro viral infections taking into account spatial characteristics of virus spreading within a culture well. The aim of the model is to understand the key mechanisms of SARS-CoV infection dynamics during the first 24hours post infection. We interrogate the model using a Latin Hypercube sensitivity analysis to identify which mechanisms are critical to the observed infection of host cells and the release of measured virus particles. In the past ten years there has been a growing interest in modeling viral infections for the study and characterization of host infection dynamics. Early mathematical models were typically based on ordinary differential equations (ODEs) and focused on extracting key parameters of the infections dynamics [1, 2] . Later on, interest moved to investigating how spatial relations affected the system dynamics, thus moving to cellular automata (CA) models [3, 4, 5] . More recently the attention has shifted to viral respiratory diseases due to the increasing danger of pandemics such as the "Spanish flu" that caused over 50 million deaths worldwide in 1918. Near pandemics such as severe acute respiratory syndrome (SARS), resulting in 8,096 infected cases and a 9.6% death rate, have provided unfortunate reminders their continued health risk and significant economic impact [6] . New research projects are being funded to investigate the different scales of threat, from the epidemiological spread of a virus in a population to its early phases of the infection in a cell culture. Among the papers published on this topic, in particular on the influenza virus [7, 8, 9] , there are Bocharov's ODE model [10] and Beauchemin's CA model [3, 4, 11] . Both models Open access under CC BY-NC-ND license. described the dynamics of the influenza infection in the upper epithelial areas of the lungs and were validated by clinical data. Although they take into account the effect of the immune system on the dynamics of the disease, describing the infection until its final outcome (more than a week after the first infection), none of them analyzed which mechanisms where critical in the dynamics of the first and second round of infection. To understand the key mechanisms of viral respiratory diseases, we focus instead on early phases of viral infection by using in vitro experiments observing lung cells up to 24 hours post infection (PI). The in vitro experiments provide measurements of virus titer, spatial characteristics of cell growth and, through green fluorescent protein (GFP) imaging, of the infection spread. Our computational model focuses on simulating the early stages of a viral infection in a population of cells plated on a culture well. The choice of a CA model was natural since the in vitro infections being studied use host lung cancer cell lines that form a fixed mono-layer in which spatially dependent aspects of infection may be present [12, 13] . We developed this computational model using the Multi-Agent System Visualization (MASyV) platform [3] . In contrast to previous models, we explicitly focus on the dynamics of virus spread on a population of cells, supported by experimental data from an in vitro model system. We also explicitly model the infectious viral particles as discrete entities, whereas in previous models the infection of cells followed simple CA rules depending on the states of neighboring cells. These viral particles are released by infected cells according to a specific function based on time post infection, and move over the well with a random walk algorithm. This representation allows us to model the mechanisms of virus spread in an environment where the virus is not confined and can also infect cells not adjacent to the infected ones. In Section 2 we describe the model design and its main features. We also describe the SARS infection experimental data used to parameterize the model. In Section 3 we present a sensitivity analysis that identifies the critical mechanisms characterizing the early phases of the infection. We also show that the model can explain the experimentally observed virus titer data and allows a deeper understanding of the infection dynamics in the in vitro experiments. The computational model is built using Beauchemin's MASyV platform. The software consists of a server providing I/O and supervisory services to the various client modules where the simulation is actually coded. Our choice to use MASyV was partially driven by flexible and powerful graphical visualization routines that facilitate comparison to images provided by the experimental collaborators. MASyV has a C-based API and is open source allowing finalized custom models to be easily shared. We discuss novel contributions and differences from the previous modules. The original module details are presented in Beauchemin et al. [3] . Our model reproduces a viral infection on a population of cells plated on a culture well. In our client we consider, as the target of the viral infection, Calu-3 cells that are a human airway epithelial cell line derived from human lung cancer. We model these host cells using a 130x130-site CA model where each site represents either one calu-3 cell or an empty space. At the beginning of the simulation each lattice site is initialized and labelled with "uninfected" or "empty" states as described below in Section 2.2. Uninfected cells are initially stochastically infected with virus through a first round of infection at the beginning of the simulation, described in Section 2.3, and once infected progress through the following states: Containing: initial infection state representing viral entry and hijacking of host cell mechanisms necessary for viral replication. Expressing: cell is actively producing and assembling virus capsids and genomes internally, but has not begun releasing virion. Infectious: Assembled virion is being released from the host cell according to the release function (Section 2.4) By examining the experimental viral titer data shown in Figure 1 we derived temporal delay of the state transition between Containing and Infectious. The ~ log 10 3 viral titer measurements at time points 0, 4 and 7 hours post infection are residual virus left over from the initial infection after washing. We set the delay for release of new infectious viral particles to 7 hours to represent the jump in observed virus titer between 7 and 12 hours. Cells release virus particles according to a viral release function, described below. Virus particles are entities that diffuse from one lattice point to another according to a random walk algorithm (Section 2.4). These free-floating virus particles stochastically infect uninfected cells in a second round of infection. The CA lattice is therefore like a tissue of immobile cells with infective virions moving over it. We represent only infectious viral particles to relate their count in the model with the viral titer measured in terms of Plaque Forming Unit per mL (PFU/mL). A PFU is a measure of the number of viral particles capable of forming plaques per unit of volume, thus only infectious particles are counted in the viral titer. The CA lattice is updated synchronously and the boundary conditions for both the epithelial and viral particles are periodic, i.e. epithelial cells and viral particles can grow or move outside the boundary of the virtual well. We use a Honeycomb neighbourhood: each lattice site is considered adjacent to six sites and cells can replicate only if a neighbouring site is empty or if it contains a dead cell. Only uninfected cells are allowed to replicate. When cells on the lattice are initialized all uninfected cells are stochastically assigned a TIME_TO_DEATH value between 1 and an arbitrary CELL_LIFESPAN variable. CELL_LIFESPAN was chosen to be large enough such that 5%, a minority of the cells, die naturally within the simulated 24 hours. CELL_LIFESPAN is not affected by infection as in vitro experimental observations showed no increase in cell death due to infection within 24 hours. In the laboratory before an in vitro infection is performed, 10 6 cells are first plated in each of the 6 culture wells present in a dish, and given time to settle down on the well surface to create a cell mono-layer, before the virus is applied to the cell culture. In the plating process, a suspension of cells is plated onto the surface of the well that adhere to the plastic in small clumps or as single cells forming "islands" that are homogeneously distributed throughout the well. Cells located on the inside of these islands (surrounded by other cells) do not replicate due to contact inhibition whereas cells on the outer edges of each island replicate until they are completely surrounded by other cells. As a result each island continues to grow until there is no adjacent empty space on the well surface. To represent this behaviour we used a set of images ( Figure 2 ) that shows how cells are distributed on the culture well just before the infection is performed. The number of cells forming each isolated island was counted for each image obtaining a distribution of island sizes with an average of 116 cells per island and a standard deviation of 74 cells per island. This data was used to generate randomly positioned clusters of cells on the grid, each with a size extracted from the measured distribution. Islands are placed iteratively until the simulated confluency matches the experimentally measured confluency. A cell line doubling time is defined as the time required for the population of cells to double in number and we used a doubling time of 48 hours for Calu-3 cells (Ardizzoni Lab, personal communication). Cell density directly affects local cell growth since cells completely surrounded by other cells are not able to replicate. For this reason in our model the initial spatial distribution of uninfected cells affects the overall cell growth on the virtual culture well. In order to obtain a correct cell growth in our model we used a parameter called DIVISION_TIME that measures the time necessary to a uninfected cell to duplicate. We adjusted DIVISION_TIME to match the simulated doubling time with the experimentally-derived one. This was accomplished by analyzing cell growth of uninfected cells over 48 hours as a function of time for different values of the parameter. This led us to fix the value of DIVISION_TIME to 10 hours. Cells in the culture well are initially infected using a Multiplicity of Infection (MOI) of 2 (see section 2.7.4). The MOI is calculated as the ratio of infectious viral particles to the number of Calu-3 target cells. By definition the proportion of infected cells is given by the Poisson distribution that describes the infection process [14] . Using the experimental MOI would theoretically give us the exact percentage of infected cells at the beginning of the experiment but even though the number of particles can be measured with good accuracy not all the viral particles used in this first infection process are infective and the estimate is derived using VeroE6 cells. This leads to the definition of an effective MOI, that is the MOI given by the number of particles able to truly infect the Calu-3 host cells. Initially viral particles infect a proportion of the plated Calu-3 cells before being washed away. We use the resulting proportion of infected cells estimated through the standard MOI definition as a starting point of our simulation. After cells are initialized on the lattice, they are assigned an infected state according to Equation 2, which describes the probability of at least one viral particle entering a cell. In the equations below, n represents the number of virus particles and MOI is the multiplicity of infection used in the experiment. Equation 2, the probability of a cell being infected by 1 or more virus particles, is derived from Equation 1, the probability of being infected by a single virus particle, commonly used to describe viral infection as a Poisson process [14] . Although the expected MOI is experimentally known, as mentioned in the "Experimental Data" section it is often over-estimated. We add the free parameter INFECTIOUSNESS for two purposes: 1) to scale MOI over-estimation and 2) predict the number of initially infected cells when this data is not available. In the sensitivity analysis we discuss the importance of this parameter. MOI n * e MOI n! (1) After cells transition to an infectious state virus particles are released with a probability described by a sigmoidal function shown in Figure 4 . The probability that one or more infectious viral particles are released in each time step is given by V(t) (Equation 3 ). In this equation t represents time since the cell was infected and parameters A, B, and C are derived by fitting an experimentally-derived viral release curve produced by Ka-Wai et al. [15] . Since we had no direct measure of the amount of virions per cell released by SARS infected cells we rescaled the sigmoidal function with the free parameter VIR_RELEASE. When V(t) > 1 an infectious cell will release one or more virions with a probability of 1. Decimal digits less than one are treated as a probability for an additional virus particle to be released. Infectious virions diffuse in the virtual well according to a simple random walk: for each virion in the virtual well at each time step of the simulation we perform a number of Random Walk steps given by the free parameter NUM_DIFF_STEPS. Each viral particle has an equal chance to move to one of the six neighboring lattice sites and at each time step performs NUM_DIFF_STEPS movements on the lattice. The need of using a free parameter for the diffusion of virions was given by the lack of data regarding the virion diffusion coefficient on the well with specific conditions of the media used to cover the cell culture. Viral titer measures only infectious particles present in the media using Plaque Forming Units (PFU) per ml. For this reason we model only infectious particles [14] . In addition it is more computationally efficient to release and diffuse only the infectious viral particles. Once released, virus particles may infect cells with an uninfected state located at the same lattice site. We assume each viral particle independently infects a nearby cell with a probability following the binomial distribution. For each lattice site the probability P second round (Equation 4) is used to calculate whether or not a given uninfected cell is infected by n viral particles given N total local virus particles. p BP represents the probability of a virus-receptor binding event leading to a cell's infection by a single viral particle during a given model time step. We refer to the free parameter, p BP , as BINDING_PROB below. Once a cell is successfully infected, n viral particles are removed from the local virtual media located at the cell's lattice site. For these experiments we used a virus derived from our severe acute respiratory syndrome coronavirus (SARS-CoV) wild type infectious clone in which we engineered the green fluorescent protein (GFP) in place of open reading frame 7a/b. SARS-CoV GFP stock titers were calculated using standard plaque tittering methods. Briefly, confluent monolayers of VeroE6 cells in 6 well plates were infected with serial 1:10 dilutions (usually 10e-1 to 10e-6) of stock virus for 1 hour at 37 o C. The monolayers were covered with a solution of 0.8% low melting point agarose (Seachem), 1X minimal essential media high glucose (MEM, Invitrogen), 10% fetal clone II (Hyclone), and 1X antibiotic/antimycotic (Invitrogen), which solidifies trapping the virus to allow cell to cell viral spread but not release of virions into the media. Plates were incubated at 37 o C for 36 hours, stained with neutral red for 2 to 5 hours, the stain removed and the plaques (holes in the monolayer generated when viruses kill the host cells) are visualized and counted on a light box. Stock titers were calculated as plaque forming units (pfu) per mL. Calu-3 cells were plated at a density of 1x10 6 per well in 6 well plates in MEM containing 10% fetal clone II and 1X antibiotic/antimycotic. Cells were incubated at 37 o C at 5% CO 2 levels for 2 days prior to infection with a media change 24 hours post plating. SARS-CoV GFP stocks (7.5x10 7 ) were used to infect each well at a multiplicity of infection of 2 or the addition of 2 infectious virions per cell in each well. At each time point, 100uL of media was harvested from each well to titer via plaque assay. (Number of cells per well assumed to be 2x10 6 at 2 days post plating.) Images shown in Figure 2 were taken using phase contrast microscopy images were taken at standard exposure times at 10X magnification. At the indicated times post infection, images of the infected and mock-infected Calu-3 cells were taken in living cells real time using a fluorescent light to excite GFP [16] At each time point, three to five fields of SARS-GFP infected cells were assessed using ImageJ (http://rsbweb.nih.gov/ij/). In ImageJ, we gated the light signal to generate spots in each cell and used the spot counting algorithm determine the total number of cells per field. GFP positive cells were then counted and averaged. For each candidate parameter set, θ, a simulation fitness (Eq. 5) was calculated based on a comparison of two experimental measurements: virus titer and proportion infected cells. Both components of simulation fitness, F(θ), are normalized by a maximum error (ME) term to balance their contribution. F titer = log 10 VT exp(t) log 10 VTsim(t) 2 (6) Randomized weight bootstrapping (1000 iterations) was used to determine ME titer and the derivation of ME GFP is described below. Equation 6 shows F titer is the sum of squares difference between experimentally derived virus titer, VT exp , (Figure 1 ) and scaled virus titer produced by the simulation, VT sim . Only the last three time points (7, 12, 24 hours) are compared because earlier time point titer values were skewed by residual virus. Simulation output, VT sim , is scaled by α to account for differences in the number of simulated cells versus the number of in vitro cells. Using a lattice size of 130x130 an average of 7,636 total starting simulated cells were produced. The virus titer was scaled by a factor of 10 6 /7,636 because the experimental virus titer was produced from a population of 10 6 in vitro cells. As described in the Experimental Data Section 2.7, in vitro GFP measurements were used to quantify the proportion of infected cells at 12 hours post infection. This additional biological information was used to constrain our model's parameter space. The F GFP portion of simulation fitness (Equation 7) is the absolute value of the difference between simulated proportion of infected cells I sim to experimentally measured proportion infected, I exp . Since I exp was measured to be .11 (on average 11 of 100 cells are infected) at 12 hours, a maximum error, ME GFP , of .89 was used to normalize F GFP between zero and one. To assess the importance of model free parameters to simulation outcome, a Latin Hypercube sampling (LHS) sensitivity analysis was performed [17, 18] . This two-tiered sampling method was implemented by dividing each free parameter's range into 1000 intervals. Intervals for each parameter were shuffled without replacement and then randomly sampled using a uniform distribution resulting in 1000 samplings of the parameter space. All possible probability parameter values were considered for INFECTIOUSNESS AND BINDING_PROB between 0 and 1. Maximum values for NUM_DIFF_STEPS and VIR_RELEASE were arbitrarily chosen (12 and 30 respectively). A Spearman rank correlation was then used to measure statistical dependence, Rho, between free parameter and single run simulation fitness. The following are free parameters (described above) contained in our model and assessed using LHS: We also examined how free parameters affect simulation output with time. Using LHS sampling, we tested relations between free parameters and two simulation outputs: virus titer and % of infected cells. Figure 6 (A) shows that INFECTIOUSNESS and VIR_RELEASE are significant over 24 hours of simulated infection and that their correlations start to diverge at ~12 hours. Again, BINDING_PROB and NUM_DIFF_STEPS do not show a significant correlation with virus titer. However, parameters necessary for the second round of infection do impact the proportion of infected cells shown in Figure 6 (B) where we see an increase in BINDING_PROB and VIR_RELEASE correlation between 7 and 12 hours, corresponding to a decrease in the INFECTIOUSNESS correlation. This decrease can be attributed to the second round of infection in which a greater population of cells has been infected. The simultaneous increase in correlation of BINDING_PROB and VIR_RELEASE indicates that additional rounds of infection are playing a significant role in viral spread. Although simplistic compared to in vivo model systems, the interpretation of in vitro experiments is still confounded by biological complexity and disparate data types. Explanatory models are critical for understanding and hypothesis generation. This agent-based modeling framework may be used to investigate first and second round of infection mechanisms using free parameters able to be tuned to allow the model to incorporate disparate types of experimental data. We also take into account spatial aspects of infection, including biases in culture well cell growth and diffusion of infectious virions. Virus titer data and GFP infectivity data from a SARS infection of Calu-3 cells is used as an example to illustrate the model's capacity to interpret experimentally derived data. LHS sensitivity analysis indicates that a small population of cells is initially infected and that additional rounds of infection are responsible for virus titer measurements. A significant relationship between INFECTIOUSNESS and both the simulation fitness and simulation outputs (virus titer and proportion infected cells) indicates the importance of this parameter on the resulting infection dynamics. This result also demonstrates the importance of the accurate cell and infectious virus particle counts for the MOI calculation. Finally our model highlights the importance of intracellular processes leading to virus release. One possible future step is to include additional detail regarding intracellular processes of virus replication and move to a multi-scale spatio-temporal model. Future work is planned to incorporate microarray data and make predictions regarding host response and expression in order to examine connections between infection state and signaling in the immune response. Simulated annealing has been used to identify free parameters that fit the described model to virus titer data and may be used to predict the number of infected cells or other un-measured data types to support experimental modeling efforts. These off-line predictions could then be used to interpolate a single-cell function representing host response postinfection. We also plan to train the model with data from multiple virus strains to investigate how virus population and host response dynamics differ. Finally, we also plan to investigate the effects of initial spatial distribution of infected cells on viral pathogenesis for multiple virus strains. All authors must sign the Transfer of Copyright agreement before the article can be published. This transfer agreement enables Elsevier to protect the copyrighted material for the authors, but does not relinquish the authors' proprietary rights. The copyright transfer covers the exclusive rights to reproduce and distribute the article, including reprints, photographic reproductions, microfilm or any other reproductions of similar nature and translations. Authors are responsible for obtaining from the copyright holder permission to reproduce any figures for which copyright exists. HIV-1 dynamics in vivo: virion clearance rate, infected cell lifespan, and viral generation time Modelling viral and immune system dynamics A simple cellular automaton model for influenza A viral infections Probing the effects of the well-mixed assumption on viral infection dynamics Structured Model of Influenza virus replication in MDCK cells SARS cases with onset of illness from 1 Kinetics of Influenza A Virus Infection in Humans Antiviral resistance and the control of pandemic influenza: The roles of stochasticity, evolution and model details Modeling the Viral Dynamics of Influenza A Virus Infection Mathematical model of antiviral immune response III. Influenza A virus infection Modeling Influenza Viral Dynamics in Tissue. Artificial Immune Systems Multi-scale modelling in computational biomedicine Modeling Dynamic Systems with Cellular Automata, Chapter 21 Role of ATP in influenza virus budding Severe acute respiratory syndrome coronavirus infection of human ciliated airway epithelia: role of ciliated cells in viral spread in the conducting airways of the lungs A methodology for performing global uncertainty and sensitivity analysis in systems biology Critique of and limitations on the use of expert judgements in accident consequence uncertainty analysis This work was made possible by funding from the National Institute of Allergy and Infectious Diseases, NIH, Department of Health and Human Services contract # HHSN272200800060C Thanks to Casey Long for his work in plating Calu-3 cells for SARS-CoV infection experiments. Also thanks to Dr.Andrea Cavazzoni of the Unit of Experimental Oncology (University of Parma) for the useful discussions.