Deterministic model

Based on the model in [24], we propose a new model which includes the virus and DIP productions. As shown in Fig 1, in the model, we consider free natural infectious virus, denoted by , and defective interfering particles (DIPs), denoted by where is a vector which represents a spatial location in a two-dimensional domain [0, x_{1 max}] × [0, x_{2 max}]. Also, there are six types of cells: uninfected cells, cells infected only by natural viruses but not in the period of virus production, cells infected only by natural viruses and in the period of virus production, cells infected by DIPs only, cells infected by DIPs and natural viruses but not in the period of virus production, and cells infected by both DIPs and natural viruses as well as in the period of virus production. The numbers of the respective cells are denoted by C, C_V, , C_D, C_VD and , respectively.

Download:

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

Fig 1. Schematic diagram of the virus-DIPs system.

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

Considering infection by DIPs of cells late in the replication cycle is too late to affect the production of the virus, we assume that DIPs cannot infect the cells [24]. There are two age categories for each of the C_V and the C_VD cells in our model. Ultimately, the DIPs are produced only by the mature cells, and virus particles are produced by both and cells. While C_V cells can get to maturity by themselves and produce virus, C_D cells cannot produce DIP unless they are co-infected by the virus and become C_VD and get to maturity. The latter models the situation that DIPs cannot replicate unless they co-infect a cell with a wild-type virus.

The following two equations are for modeling the dynamics of the virus and DIP: (1) where ∇² is the Laplacian operator, describing the virus and DIP diffusion.

We assume that cells are not moving in the spatial domain and we can model the dynamics of the cell densities by the following system: (2) where is the total density of all cells.

Our model Eqs (1) and (2) is different from that of [24] in several ways: (i) there are two age categories for C_VD cells in our model, but there is no age structure for C_VD cells in [24]; (ii) the mature cells can produce virus in our model, but the C_VD cells in [24] cannot produce virus; (iii) C_D cells cannot recover to be uninfected cells in our model, but they can recover in [24]; (iv) our parameters γ₁ and γ₂ can be different, but they are the same in [24].

Stochastic models in two different scenarios of virus production

Since the outcome of the model with DIPs is sensitive to the competition between viruses and DIPs, different kinds of perturbation to the production of viruses and DIPs may contribute to a huge change in the probability distribution of the outcome. The study in [29] suggested that there are two scenarios of virus production, which can create different kinds of perturbations to virus production:

1. Scenario 1: infected cells produce virus and DIPs through cell bursting;
2. Scenario 2: infected cells keep producing viruses and DIPs continuously.

However, these two scenarios cannot be distinguished by our deterministic PDE model [29] as both models with different scenarios have identical mean-field kinetics. In this study, we built a stochastic model and developed an efficient simulation method to examine the effects on the spatial distribution of viruses under different scenarios.

Due to the high computational cost of the spatial stochastic model, there are not many studies considering the effects of different scenarios for virus production on the spreading speed and distribution of the virus. To improve the computational efficiency, here we simulate our model with Spatial Stochastic Simulation Algorithm (SSA) [32], which is a method to generate an exact sample from the probability mass function that is the solution of the chemical master equation.

In SSA, we consider the spatial domain as a two-dimensional square with length L. The domain is partitioned into N_c × N_c identical compartments that are uniform squares with length h = L/N_c. The subsystem in each compartment is assumed to be homogeneous. The same types of particles and cells in different compartments are treated as different species; for example, we denote by V_i,j the virus level in the compartment at location (i, j) and consider . Diffusion is treated as a reaction in which a molecule jumps to one of its neighboring compartments at a constant rate. Then with no-flux boundary conditions (or other conditions which depend on the experimental setting), diffusive jumps obey the following chain reactions for each j ∈ {1, 2, ⋯, N_c}: where ρ₁ = d_V/h². We assume that D_i,j has similar chain reactions with ρ₂ = d_D/h². We define the propensity function for the jumps, for example, at the location (i, j), for the four types of jumps (L: left, R: right, U: up, D: down) of virus: , , , and . At the boundary, some jumping directions will not be considered for no-flux boundary conditions. For reactions, we assume that only molecules in the same compartment can react with each other.

Different scenarios of virus production will contain different sets of reactions. In the first scenario, the reactions in the (i, j) compartment are as follows:

In the second scenario, the reactions (the first three rows of the previous scenario) are the same as the first one except for the production of viruses and DIPs. That is, we replace the last row by the following:

The domain and multiple interfaces for different reactions

Consider a general reaction-diffusion system of S species and M chemical reactions and diffusion in 4 directions in a two-dimensional domain Ω, which is partitioned into N_c regular compartments of width h. Let N_s(k, t) represent the amount of the s-th species in the k-th compartment at time t. Each compartment is small enough so individuals in it can be assumed well mixed.

The subdomain in which we employ the compartment-based regime for the j-th reaction at time t is denoted by , and the other part of Ω that employs PDE is represented by . contains all compartments in which the amount of at least one of the reactants in the j-th reaction is below the threshold value Θ. To be specific, assume that reactants of the j-th reaction are {S₁, S₂, ⋯, S_m}. If (3) then the k-th compartment is assigned to the stochastic domain , otherwise to the PDE domain . The parameter Θ is introduced due to the inherent cut-off in experiments, such as the technical minimum detectable level of RFP/GFP. We tested different values of Θ and selected the most appropriate one. In our algorithm, interfaces are adaptive. Domain division and multiple interfaces are updated every Δt_I. Fig 2 shows a one-dimensional illustration of the approach stated above.

Download:

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

Fig 2. An illustration of the domain division and interface of the j-th reaction.

Here we show an example with two reactants. The domain is modeled by PDE and the domain is modeled by compartment-based SSA. The amount of each species in a compartment in is . If any reactant amount is below the threshold Θ, then that compartment is part of . Individuals can move between the boundary compartment of and the pseudo-compartment in . In the two-dimensional case, diffusion takes four directions: up, down, left and right.

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

It’s worth noting that and can be the same if there is an inclusion relationship between the sets of reactants in the j₁-th and j₂-th reaction. In fact, the number of non-coincident interfaces is no more than the total number of species in the system. Therefore, compared with a single interface, multiple interfaces can capture stochastic fluctuations more accurately without increasing too much computation costs.

The pseudo-compartment method

In this section, we will outline the pseudo-compartment method [35], which is the basis of our algorithm. In [35], Yates et al. introduced the pseudo-compartment method for diffusion. On this basis, we propose the possibility of multiple adaptive interfaces.

Consider a reaction-diffusion system of S species, M chemical reactions and diffusion in the four directions of the cross in a 2D domain Ω. In our algorithm, the PDE region varies for each reaction. So instead of just dividing the PDE based domain, we discretize the whole domain, Ω, into a regular grid with spacing Δx. We consider the density of each species. For the j-th reaction at time t, the PDE numerical solution is updated for all grid points lie in . Diffusion terms are treated in a similar way, but employing the implicit Euler method. A zero-flux boundary condition is implemented in , including domain boundaries and interfaces. Flux at the interface is implemented in the compartment-based regime.

The compartment-based regime evolves from the Gillespie algorithm (SSA). Consider the propensity function of reactions and diffusion, α_i,j(t), for compartment . α_i,j(t)dt represents the probability that the j-th reaction (for j ∈ {1, ⋯, M + 4}, including diffusion) occurs in C_i during the small time interval [t, t + dt].

The coupling is implemented with a pseudo-compartment, C₋₁, presented for diffusion between the deterministic and stochastic domains. This is a compartment adjacent to the interface but within deterministic domain , where j ∈ {M + 1, ⋯, M + 4}, representing diffusion (four directions of the cross). In order to correctly model the flux over the interface, individuals in the boundary compartment in can jump into the pseudo-compartment with the usual diffusive rate, and vice versa. The amount of each species within the pseudo-compartment is calculated through direct integration of the PDE, (4) where p_s(x, t) is the PDE solution of density of the s-th species. Then the propensity function for jumping from the pseudo-compartment to the adjoining compartment in is given by (5)

The Gillespie’s direct method [36] is used to simulate the time evolution of stochastic regime. The time interval for next reaction, τ, is determined by: (6) where r₁ is a random variable uniformly distributed in (0, 1). We then use the second random number r₂ to find the corresponding reaction or jump. The algorithm then checks whether the closer update time is for PDE or SSA. If t + τ < t_P, then the update is for SSA and t = t + τ; otherwise it is for PDE and t = t_P, t_P = t + Δt_P.

Moving interface

The multiple interfaces are updated with time step Δt_I, by recomparing amounts in a compartment of all reactants of each reaction with the threshold Θ. Similar to [37], after the interfaces are updated, we need to keep numbers in the stochastic domain as integer values, but we cannot simply get rid of the fractional parts. Suppose the compartment C_k is moved from the PDE domain to the stochastic domain, and the fractional part is (7) where {⋅} is the fractional part function. P is used as the probability that an additional individual is kept in this compartment. We then take a uniform random number r ∈ [0, 1]. If r < P then we place the individual in compartment C_k; otherwise it is placed in the deterministic domain.

The pseudocode for our algorithm is given in Algorithm 1.

Algorithm 1 The hybrid algorithm for reaction-diffusion systems.

1. Initialize the time, t = t₀ and set the final time, T. Specify the PDE-update time step Δt_P and initialize the next PDE time step to be t_P = t + Δt_P. Specify the interface-update time step Δt_I and initialize the next interface-update time step to be t_I = t + Δt_I.

2. Specify the PDE spacial step Δx and the compartment width h. Initialize the amount of each species in each compartment, N_s(k, t) for k ∈ {1, …, K} and specify the threshold Θ. Compute the density, p_s(x, t) = N_s(k, t)/h for PDE grid points.

3. Determine the initial interface for each reaction j, j ∈ {1, 2, ⋯, M}:

 (a) Find all k such that , where S_j contains all species involved in reaction j, then the k-th compartment is part of the stochastic domain , and otherwise part of the PDE domain .

 (b) All compartments adjacent to (no diagonal angles) are regarded as pseudo compartments.

4. Determine the time for the next ‘compartment-based’ event according to the Gillespie algorithm, t_C = t + τ.

5. If min{t_C, t_P, t_I} = t_C then the next compartment-based event occurs:

 (a) Determine which event occurs according to the Gillespie algorithm.

 (b) If the event is moving from stochastic domain to a pseudo compartment, C₋₁, then for the corresponding (s, k), N_s(k, t + τ) = N_s(k, t) − 1 and . Here, is an indicator function that takes the value 1 when x ∈ A and 0 otherwise.

 (c) If the event is moving from a pseudo compartment C₋₁ to stochastic domain and p(x, t) > 1/h for all x ∈ C₋₁, then N_s(k, t + τ) = N_s(k, t) + 1 and .

 (d) Update the density for the pseudo compartment.

 (e) Update the current time, t = t_C.

6. If min{t_C, t_P, t_I} = t_P then the PDE domain is updated:

 (a) Apply backward Euler for diffusion terms and forward Euler for reaction terms.

 (b) Update the density for the pseudo compartment.

 (c) Update the current time, t = t_P and set t_P = t + Δt_P.

7. If min{t_C, t_P, t_I} = t_I then the interfaces are updated, similar to step 3:

 (a) For each reaction, find all k such that , where S_j contains all species involved in reaction j, then the k-th compartment is part of the stochastic domain , and otherwise part of the PDE domain .

 (b) All compartments adjacent to (no diagonal angles) are regarded as pseudo compartments.

 (c) For the compartment C_k that change from PDE domain to stochastic domain, let . Take a random number r_s ∈ [0, 1].

  • If r_s < P_s then N_s(k, t_I) = ceil (N_s(k, t)) and p_s(x, t_I) = p_s(x, t) − (1 − P_s)/ Area for ;

  • otherwise, N_s(k, t_I) = floor (N_s(k, t)) and p_s(x, t_I) = p_s(x, t) + P_s/ Area for .

 (d) Update the current time, t = t_I and set t_I = t + Δt_I.

8. If t ≤ T, return to step 4.

 Else end.

Parameter selection

We consider the spatial domain as a two-dimensional square with length L = 2.552mm, which is the same as the experimental data; for the PDE numerical scheme, we apply the central difference scheme to discretize the Laplace operation with Δx = Δy = 0.058mm; for the temporal numerical scheme, we use the backward Euler method for the Laplace operation and forward Euler method for the other terms with time step Δt = 0.01h. In the SSA approximation, the domain is partitioned into square compartments with dimension h × h = Δx × Δy.

Diffusion coefficients of virus and DIP are set to be 2.38 × 10⁻⁶cm² in [31] while the decay rate is 4.0 × 10⁻⁵s⁻¹. As the diffusion rate varies according to the environment and plays a vital role in spatial distribution, we increased the former d_V = d_D = 2.38 × 10⁻³mm²/h to provide qualitative agreement with experimental data and left the latter unchanged δ_V = δ_D = 0.144h⁻¹.

In [27], the rate of virus production is expressed as the product of the number of viruses released per cell after packaging and the rate at which each cell produces viruses. Therefore α₁ = 758.045 × (68.503 × 10^±2d⁻¹) = 2163.682 × 10^±2h⁻¹, and α₃ = 38.259 × (21.782 × 10^±2d⁻¹) = 34.723 × 10^±2h⁻¹. The wide range of parameters allows us to choose a suitable value to obtain a good agreement with the experimental results qualitatively. So we set α₁ = 6.491h⁻¹ and α₃ = 69.446h⁻¹. Since DIPs may exhibit a replication advantage over infectious viruses [3], we assumed α₂ = α₃/10 in this work.

Same as [27], the intrinsic rate of uninfected cell proliferation is α_C = 15.217d⁻¹ = 0.634h⁻¹. But the cellular carrying capacity of proliferation varies depending on the experimental environment. We let K = 3.505 × 10⁵ × h² cells/compartment to achieve qualitative agreement with experimental data, where h² is the compartment area.

The rate of maturation of C_V cells into cells is 9.863 × 10^±2 d⁻¹ in [27]. We slightly increase it to ν₁ = ν₂ = 0.205h⁻¹ because mature infected cells are observed later in experiments. β₁ and β₂ are considered as the death rate of and that of respectively, which are 2.426 × 10^±2d⁻¹ in [27]. We take β₁ = β₂ = 0.05h⁻¹ in simulations.

Virus and DIP infection rate is 2.45 × 10⁻¹⁰ = 1.02 × 10d⁻¹ = 1.02 × 10⁻¹¹h⁻¹ in [27], which is relatively small. Different experiments and higher cell density may lead to a larger infection rate. Hence we set γ₁ = γ₂ = 4 × 10⁻⁴h⁻¹.

The infected cell death rate is 5.91 × 10⁻²h⁻¹ in [31], which is used as death rates for all cells in our simulations.

All parameters are listed in Table 1. It is worth noting that our set of parameters can guarantee that species in the system without DIPs will coexist in the following simulations. A detailed proof is provided in S1 Appendix.

Download:

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

Table 1. Parameters used in the simulations.

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

Interpretation of experimental data

To validate the usefulness of our model, we compare the simulation results with the experiment data published in [3], the latter is composed of time series of images obtained via Microscopy from the co-propagation of infectious and defective viruses in a population of biological cells. These co-infection experiments were initiated with the same virus inputs (MOI 30) but different DIP inputs (namely MOI 0,1,10 and 84). and microscopy images were taken at 7 hours, 13 hours, 19 hours and 25 hours post infection. The DIP expresses a green fluorescent protein (GFP) and the wild-type virus expresses a red fluorescent protein (RFP). There are three to five time series for each of the RFP intensity and the GFP intensity. Each image has size of (2200, 2200) with diameter of 1.16 μm pixel. The scale bar is 0.5 mm.

Fluorescent protein labeling is usually used for qualitative purposes, and there is no linear relationship between brightness and intensity. Therefore the experimental images only provide a reference for virus expression in simulations.

Since in the following simulations we employed the compartment-based method while experiments provide scatter diagrams, we have done some preprocessing to compare them with the computer simulation results. Fig 3A is a representative experiment figure. We extracted the red single channel (the virus is expressed) and filtered noise, as shown in Fig 3B. We then did morphological transformations (dilation followed by erosion) to close small holes inside the objects. Therefore Fig 3C maintains the critical features of virus expression in experiments and is more approximate to compartment-based.

Download:

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

Fig 3. A representative example of interpreting the preprocessing of experimental figures.

A: An original experiment figure. B: The red single channel is extracted (virus expression) and denoised. C: Morphological transformations (dilation followed by erosion).

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

Dynamics and pattern formation of virus expression

We first study the PDE model in Eqs (1) and (2). Fig 4 shows the spatial distribution of and in a 2D domain at time t = 25h with different initial DIPs, as well as images of experimental data at the same time in [3]. When there is no DIP, viruses are uniformly radially distributed in both simulations and experiments. When the initial condition is C_VD(0) = 200 for the PDE model, then viruses are distributed in a ring while the DIPs are radially distributed in the center. Compared with the experimental results under similar conditions, the PDE model shows good agreement in the speed of infection with experimental data but the spatial distribution does not match and the patchiness of is not observed in PDE simulations.

Download:

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

Fig 4. Spatial distributions of virus and DIP in cells in PDE simulations and experiments.

A: Spatial distributions of virus in cells () and DIP in cells () in PDE simulations and representative experiment with initial DIP equal to 0. B: Spatial distributions of virus in cells () and DIP in cells () in PDE simulations and representative experiment with initial DIP equal to 84.

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

In stochastic simulations, the same types of particles and cells in different compartments are treated as different species, for example, for V, denoted by . The initial condition is and C_i,j(0) = 1000 for all (i, j), C_Vi,j(0) = C_VDi,j(0) = 0 for all (i, j) except the midpoint C_V22,22(0) = 100, and C_VD22,22(0) varies from 0 to 400.

Two scenarios are considered in simulations and compared with experimental results:

1. Scenario 1: infected cells produce virus and DIPs through cell bursting;
2. Scenario 2: infected cells keep producing viruses and DIPs continuously.

Fig 5 shows the evolution of viruses in cells and DIPs in cells without initial DIP from t = 17h to t = 25h. Row 1 and 2 show the spatial distribution of matured infected cells and , which is proportional to viral expression and DIP—virus expression in Scenario 1 and Scenario 2, at time t = 17 h and t = 25 h respectively. The experimental observation has an inherent threshold, and images have been denoised; therefore, we also introduced a cut-off for simulation data. Namely the amount of cells is set to be zero if it is less than the cut-off value 50, which is also applied to all the following simulations. The last column is the evolution of a representative experiment without initial DIP. When there is no DIP in the system initially, there is no DIP growth, and the virus growth is radially symmetric and flat in both scenarios and experiments.

Download:

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

Fig 5. Dynamics of virus and DIP in cells in 2 scenarios simulations and representative experiment with initial DIP equal to 0.

A: spatial distribution of virus in cells () and DIP in cells () in Scenario 1 (infected cells produce virus and DIPs through cell bursting) at time t = 17h and t = 25h; B: spatial distribution of those in Scenario 2 (infected cells keep producing viruses and DIPs); C: the representative experimental results.

https://doi.org/10.1371/journal.pcbi.1011513.g005

In experiments, the radial symmetry disappears as the initial amount of DIP increases. In fact, patchy formation of cells infected by virus is sensitive to the dose of DIP. It can be observed even with a small initial dose of DIP (Fig 6). When the initial DIP is raised from 10 to 84 in experiments, the majority of the virus at the end is concentrated (see Fig 7 and S1 Fig). Simulations show similar results, but Scenario 1 shows a much higher probability of forming a pattern than S2 and matches the experimental data better. In scenario S1, infected cells produce viruses and DIPs through cell bursting, and then viruses and DIPs diffuse, which leads to a more accidental position; while in Scenario 2, infected cells keep producing viruses and DIPs, meaning the location of those cells will continuously produce virus and DIPs. Hence the spatial distribution is more centralized rather than patchy.

Download:

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

Fig 6. Dynamics of virus and DIP in cells in 2 scenarios simulations with C_VD22,22(0) = 4 and representative experiment with initial DIP equal to 1.

A: spatial distribution of virus in cells () and DIP in cells () in Scenario 1 (infected cells produce virus and DIPs through cell bursting) at time t = 17h and t = 25h; B: spatial distribution of those in Scenario 2 (infected cells keep producing viruses and DIPs); C: the representative experimental results.

https://doi.org/10.1371/journal.pcbi.1011513.g006

Download:

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

Fig 7. Dynamics of virus and DIP in cells in 2 scenarios simulations with C_VD22,22(0) = 40 and representative experiment with initial DIP equal to 10.

A: spatial distribution of virus in cells () and DIP in cells () in Scenario 1 (infected cells produce virus and DIPs through cell bursting) at time t = 17h and t = 25h; B: spatial distribution of those in Scenario 2 (infected cells keep producing viruses and DIPs); C: the representative experimental results.

https://doi.org/10.1371/journal.pcbi.1011513.g007

Spread rate of virus

To quantify the spread characteristics of viral expression under stochastic effects, we define the virus radius as: (8) where C_V(x, y, t) represents the amount of cells infected by virus at grid point (x, y) at time t. Since a certain amount of virus expression is required to be observed in the experiment and noise is filtered, we also set a cut-off θ = 50 for computing area, i.e. (9) For experimental data, a detailed illustration is in Fig 3. Fig 8 shows the radius of virus versus time 9 ≤ t ≤ 25 (h) with different initial DIP inputs in 2 scenarios simulations and experiments. We can see the radius keeps increasing and viruses keep spreading in all cases. Whereas as initial DIPs increase, in both experiments and simulations, the growth rate of radius goes down, which is due to the inhibitory effect of DIPs on viruses. On the other hand, the Scenario 1 can better match the experimental results, both in terms of the dynamics and the level of fluctuations.

Download:

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

Fig 8. Radius of virus against time with different initial DIP inputs in 2 scenarios simulations and experiments.

In each subfigure, curves in different colors are different simulations (or experiments) with the same initial conditions. The initial conditions for the same row are the same.

https://doi.org/10.1371/journal.pcbi.1011513.g008

We note that after a certain time, the plague radius grows linearly with respect to time for each fixed initial dose of DIPs, and studied the relationship between the virus radius growth rate and initial DIP dose. To get rid of the difference in units of initial conditions in simulations and experiments, we consider a dimensionless ratio . Since initial viruses remain the same, ρ is proportional to initial DIPs. Fig 9 shows the relationship between the virus radius growth rate and initial DIP dose intuitively. We used a logarithmic x-axis, so it is shifted by 0.01 to avoid troubles when ρ = 0. The growth rate was computed using the data points after t = 13 h to ensure in all cases we have close to linear growth in radius vs. time (slope of lines in Fig 8). We run 50 group simulations for each initial condition to compute the average. When there is no DIP in the system, the virus radius growth rates are the same in both scenarios; as initial DIPs increase, the growth rate drops dramatically, which means DIPs slow down the growth of virus particles.

Download:

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

Fig 9. The growth rate of virus radius against initial DIP inputs in 2 scenarios simulations and experiments.

The x-axis is a logarithmic scale and the errors are described by the 95% confidence intervals.

https://doi.org/10.1371/journal.pcbi.1011513.g009

Patchiness via q-statistic

Patchy spatial patterns are typically observed in the image data when the initial dose of DIP is large enough. To quantify the patchiness of the image data, we employ a standard spatial statistic called the q-statistic [38]. We choose this statistic because it is suitable for our data and convenient to implement. The definition of this statistic depends on our choice of strata, which is a decomposition of the image data. In our case, the entire image is divided into 30 sectors with an equal ratio of the angle to form 30 strata S = {L₁, L₂, ⋯, L₃₀} and the union of L_i is the whole plaque P. A visual illustration is as shown in S2 Fig. In our case, since experiments employed qualitative rather than quantitative methods, that is, we can see viral expression at all fluorescent points but the brightness of these points is not proportional to the intensity. So we convert all figures binary and M(i, j, t) denotes the brightness at the (i, j)-th pixel at time t for the image, which range is {0, 255}. For simulation results, we set a threshold for the binary transformation to approximate the threshold inherent in the experimental methods and offset the loss when denoising the experimental images. Specifically, when , M = 0 and the point is black in the image; when , M = 255 and the point is red (an example in S3 Fig). Then, the q-statistics is defined to be (10) where and |A| is the cardinality of set A. N, N_h, σ, σ_h denote the number of pixels in the entire image, the number of pixels in each stratum, the standard deviation of the entire image and the standard deviation of each stratum, respectively. This statistic is invariant under spatial scale and remains the same if the intensity of the image is multiplied by a factor.

A more intuitive formula for the q-statistic is as follows [38], here we omit the time dependence, meaning denote M(i, j, t) by M_i,j and q_t by q; The denominator of Eq (10) can be written as (11) Call those two terms the sum of squares within strata (SSW) and the sum of squares between strata (SSB) and note that the numerator of Eq (10) is exactly SSW, so (12)

So if q ≈ 1, that means the sum of squares within strata is relatively more minor, indicating in each stratum, the virus is concentrated and the sum of squares between strata is somewhat more significant, meaning the differences between strata are large, so the image would appear to be more patchy. If q ≈ 0, the variance in each stratum is large, and the differences between strata are minor, so the picture would appear to be not so patchy. The q-statistic provides a suitable and convenient way to quantify the patchiness that are visibly observed in the experimental data we studied in this work. By choosing different strata according to the data, like rectangular grid in a non-radially symmetric system and triangular mesh for surfaces, this method may be applicable to other studies of spatial temporal systems.

We study the behavior of q-statistic of cells infected by the virus at time t = 25 h, that is in simulations, when the initial dose of DIPs varies. To get rid of the difference in units, we consider a dimensionless ratio . Since we always keep the initial viruses constant, ρ is proportional to the initial dose of DIPs.

In Fig 10, the x-axis is a logarithmic scale so we shift it by 0.01 to the right to avoid trouble when ρ = 0 (initial DIP is zero). On the left, the blue line denotes the q-statistic of Scenario 1 (infected cells produce viruses and DIPs through cell bursting) while the red line denotes that of Scenario 2 (infected cells keep producing viruses and DIPs) at time t = 25 h. The errors are described by the 95% confidence intervals. On the right, we marked the q-statistic for each experimental image at time t = 25 h and plot the average for four groups of experiments. Both scenarios show the same trend as experiments. When there is no DIP in the beginning, the q-statistic is minor, meaning the spatial distribution is uniform. The q-statistic increases as the initial dose of DIPs increases. Taking into account the conclusions of the previous section, DIPs slow down the growth of virus particles and make them more patchy. But when the initial dose of DIP is large enough, we observe a drop of q-statistic. It may be caused by the domination of DIPs. To be specific, when the infection of DIPs is dominated, cells infected by the virus are distributed as sporadic patches and the majority of the domain is homogeneous (). The q-statistic is sensitive to the changes in DIPs. On the other hand, Scenario 2 shows a closer magnitude of q-statistic to experimental data while that of Scenario 1 is relatively higher. In fact, Scenario 1 leads to a higher level of amount of and larger size of patterns. When infected cells produce viruses and DIPs through cell bursting, their positions are more stochastic, and hence there is a larger probability of patchy formation, which also explains the wider confidence interval of Scenario 1.

Download:

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

Fig 10. The average q-statistics of against initial DIP inputs at time t = 25h in 2 scenarios simulations and experiments.

A: the blue (red) line is the q-statistic of Scenario 1 (Scenario 2), and the bands show the 95% confidence intervals; B: the q-statistic for each experimental image.

https://doi.org/10.1371/journal.pcbi.1011513.g010