Patch formation driven by stochastic effects of interaction between viruses and defective interfering particles

Defective interfering particles (DIPs) are virus-like particles that occur naturally during virus infections. These particles are defective, lacking essential genetic materials for replication, but they can interact with the wild-type virus and potentially be used as therapeutic agents. However, the effect of DIPs on infection spread is still unclear due to complicated stochastic effects and nonlinear spatial dynamics. In this work, we develop a model with a new hybrid method to study the spatial-temporal dynamics of viruses and DIPs co-infections within hosts. We present two different scenarios of virus production and compare the results from deterministic and stochastic models to demonstrate how the stochastic effect is involved in the spatial dynamics of virus transmission. We compare the spread features of the virus in simulations and experiments, including the formation and the speed of virus spread and the emergence of stochastic patchy patterns of virus distribution. Our simulations simultaneously capture observed spatial spread features in the experimental data, including the spread rate of the virus and its patchiness. The results demonstrate that DIPs can slow down the growth of virus particles and make the spread of the virus more patchy.

effects play an important role in virus spreading, but the stochastic effects in the virus and DIP transmissions are poorly understood.It inspires us to build a stochastic spatial model to study the interaction of DIPs and viruses and how the effect of DIPs leads to patchy patterns of virus expression observed in experiments.
In this paper, we develop and analyze a new mathematical model to study the spreading speed and the spatial pattern generated by the interaction of viruses and DIPs.To incorporate the random movements of the virus and the DIPs and the stochastic effect of the interactions due to finite number of particles, we developed a stochastic reaction-diffusion system for the virus and DIP co-infection and built a hybrid method for stochastic simulation.Our stochastic model enables also the study and comparison of two common scenarios of virus production.Our simulation results demonstrated that this model can regenerate simultaneously the patchy patterns and the spread rates observed in wet-lab experiments [3], which was not achieved in previous studies [21].

Modeling
Our new hybrid model is developed based on the deterministic reaction-diffusion model introduced by Frank [20], but has several differences and new features.Importantly, our model and simulation results capture spatial spread features in two-dimensions observed in experiments and overcome computational challenges in stochastic simulations in two-dimensional domains, while the results in [20] are for one-dimension.Furthermore, we introduced and compared two different scenarios of virus production in the stochastic simulations.
Below we describe firstly the deterministic part of our model which is a system of partial differential equations, and secondly our stochastic model that incorporates two different scenarios of virus production.

Deterministic model
Based on the model in [20], 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 V (t, x), and defective interfering particles (DIPs), denoted by D(t, x) where x = (x 1 , x 2 ) 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 * V , C D , C V D and C * V D , respectively.
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 C * V [20].There are two age categories for each of the C V and the C V D cells in our model.Ultimately, the DIPs are produced only by the mature C * V D cells, and virus particles are produced by both C * V and C * V D 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 V D 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: where ∇ 2 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: where is the total density of all cells.Our model Eqs (1)-Eqs (2) is different from that of [20] in several ways: (i) there are two age categories for C V D cells in our model, but there is no age structure for C V D cells in [20]; (ii) the mature C * V D cells can produce virus in our model, but the C V D cells in [20] cannot produce virus; (iii) C D cells cannot recover to be uninfected cells in our model, but they can recover in [20]; (iv) our parameters γ 1 and γ 2 can be different, but they are the same in [20].

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 [25] suggested that there are two scenarios of virus production, which can create different kinds of perturbations to virus production: Scenario 1: infected cells produce virus and DIPs through cell bursting; Scenario 2: infected cells keep producing viruses and DIPs continuously.
However, these two scenarios cannot be distinguished by our deterministic PDE model [25] 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) [28], 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, where ρ 1 = d V /h 2 .We assume that D i,j has similar chain reactions with ρ 2 = d D /h 2 .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: α LVi,j (t) = ρ 1 V i,j (t), α RVi,j (t) = ρ 1 V i,j (t), α U Vi,j (t) = ρ 1 V i,j (t), and α DVi,j (t) = ρ 1 V i,j (t).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 March 27, 2023 5/29 replace the last row by the following: A new hybrid method for stochastic simulation In general, the computational cost for a stochastic simulation of a system in two-dimensional domain is extremely high.To reduce the computational cost and maintain the accuracy, we built up a new hybrid method which combines the advantages of our previous works: method of operator splitting [29], and spatially coupled hybrid method with adaptive interface [30].In the new method, we use operator splitting to improve the efficiency and maintain the accuracy of the simulation; also, through this method with mixing stochastic and deterministic methods, we can apply the hybrid method for specific reactions while keeping others deterministic and hence consider only part of random effects to study which stochastic behavior plays an essential role in the pattern formation.
The hybrid method combines two classes of simulation methods for modeling the reaction processes at two different scales.To capture the advantages of the methods with different scales, we use the method in our previous work [30] to separate the spatial compartments into two types of regions with adaptive interfaces: 1) the regions with "large" numbers of molecules; 2) the regions with "small" numbers of molecules.A more precise criteria for determining "large" and "small" will be given in Eq (3).
To better adapt to the complex system, we separate the compartments for each operator independently.That is, only the number of molecules of the species involved in an operator is considered in the regional division of that operator.We then apply SSA to approximate the dynamics in the region (1), and apply the PDE approximation in the region (2).For coupling two regions, we will apply the pseudo-compartment method [31] with the adaptive interface method we used in [30] in which the locations of the interfaces between different approaches are changing according to the distribution of molecules.With the idea of operator splitting mentioned above [29], our method can provide a numerical framework for studying the spatial stochastic effect of virus transmission.Through this new tool, we will have an efficient method to gain a quantitative understanding about the spatial effect of DIPs in virus transmission.

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 Ω j C (t) ⊂ Ω, and the other part of Ω that employs PDE is represented by Ω j P (t).Ω j C (t) 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 then the k-th compartment is assigned to the stochastic domain Ω j C (t), otherwise to the PDE domain Ω j P (t).In our algorithm, interfaces are adaptive.Domain division and If any reactant amount is below the threshold θ, then that compartment is part of Ω j C (t).Individuals can move between the boundary compartment of Ω j C (t) and the pseudo-compartment in Ω j P (t).In the two-dimensional case, diffusion takes four directions: up, down, left and right.multiple interfaces I j = Ω j P (t) ∩ Ω j C (t) are updated every ∆t I .Fig 2 shows a one-dimensional illustration of the approach stated above.
It's worth noting that I j1 and I j2 can be the same if there is an inclusion relationship between the sets of reactants in the j 1 -th and j 2 -th reaction.In fact, the March 27, 2023 7/29 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 [31], which is the basis of our algorithm.In [31], 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 Ω j P (t).Diffusion terms are treated in a similar way, but employing the implicit Euler method.A zero-flux boundary condition is implemented in Ω j P (t), 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 C i ⊂ Ω j C (t). α i,j (t)dt represents the probability that the j-th reaction (for The coupling is implemented with a pseudo-compartment, C −1 , presented for diffusion between the deterministic and stochastic domains.This is a compartment adjacent to the interface but within deterministic domain Ω j P (t), 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 Ω j C 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, 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 Ω j C (t) is given by The Gillespie's direct method [32] is used to simulate the time evolution of stochastic regime.The time interval for next reaction, τ , is determined by: where r 1 is a random variable uniformly distributed in (0, 1).Use the SSA with the second random number r 2 to find the corresponding reaction or jump.

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 [33], after the interfaces are updated, we need to keep numbers in the stochastic domain are 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 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.

Parameter estimation
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 −6 cm 2 /h in [27] while the decay rate is 4.0 × 10 −5 s −1 .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 −3 mm 2 /h to match the experimental data and left the latter unchanged δ V = δ D = 0.144h −1 .
Same as [23], the intrinsic rate of uninfected cell proliferation α C = 15.217d−1 = 0.634h −1 .But the cellular carrying capacity of proliferation varies depending on the experimental environment.We let K = 3.505 × 10 5 × h 2 cells/compartment to match the experimental data, where h 2 is the compartment area.
Virus and DIP infection rate is 2.45 × 10 −10 d −1 = 1.02 × 10 −11 h −1 in [23], which is relatively small.Different experiments and higher cell density may lead to a larger infection rate.Hence we set The infected cell death rate is 5.91 × 10 −2 h −1 in [27], 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 Appendix.

Interpretation of experimental data
The experiment data published in [3] 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.5mm.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.

Dynamics and pattern formation of virus expression
We first study the PDE model in Eqs (1)-Eqs (2)  In stochastic simulations, the same types of particles and cells in different compartments are treated as different species, for example, for V , denoted by 0) = 0 and C i,j (0) = 1000 for all (i, j), C V i,j (0) = C V Di,j (0) = 0 for all (i, j) except the midpoint C V 22,22 (0) = 100, and C V D22,22 (0) varies from 0 to 400.
Two scenarios are considered in simulations and compared with experimental results: March 27, 2023 13/29 Scenario 1: infected cells produce virus and DIPs through cell bursting; Scenario 2: infected cells keep producing viruses and DIPs continuously.
Fig 5 shows the evolution of the virus without initial DIP.The first two rows are time series of amounts of matured infected cells C * V and C * V D , which is proportional to viral expression and DIP-virus expression in Scenario 1 from time t = 9h to 25h.The third and fourth rows are those in Scenario 2. 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 row 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.
In experiments, the radial symmetry disappears as the initial amount of DIP increases.In fact, patchy formation 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,Fig A.1). 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.

Spread rate of virus
To quantify the spread characteristics of viral expression under stochastic effects, we define the virus radius as: 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.
Area(C V (x, y, t)) = (x,y) For experimental data, a detailed illustration is in Fig 3. 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.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 ρ = C V D (0) C V (0) .Since initial viruses remain the same, ρ is proportional to initial DIPs.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 = 13h 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 for computing 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.

Patchiness via q-statistic
Patchy spatial patterns are typically observed in the image data when the initial dose of DIP is large enough.Therefore, we quantify the patchiness of image data by the q-statistic, which is a standard spatial statistic used to measure spatial stratified heterogeneity [34].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 March 27, 2023 16/29  .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 3).Then, the q-statistics is defined to be where M t A = (i,j)∈A M (i,j,t)

|A|
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 [34], here we omit the time dependence, meaning denote M (i, j, t) by M i,j and q t by q; The denominator of Eqs (10) can be written as March 27, 2023 19/29 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 Eqs ( 10) is exactly SSW, so 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.
We study the behavior of q-statistic of cells infected by the virus at time t = 25h, that is C * V (25) in simulations, when the initial dose of DIPs varies.To get rid of the difference in units, we consider a dimensionless ratio ρ = C V D (0) C V (0) .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 green line denotes that of Scenario 2 (infected cells keep producing viruses and DIPs).The 95% confidence intervals are also presented respectively.On the right, we marked the q-statistic for each experimental image at time t = 25h 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.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.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.

Discussion
DIPs can co-infect a cell with viable viruses and interfere with virus production [1,3].However, the mechanism by which DIPs affect the spatial distribution of virus expression is still unclear.
In this work, we constructed a PDE model to describe the interaction between viruses and DIPs in a two-dimensional domain.Moreover, to study the stochastic effect on spatial dynamics of the virus spreading and patchy formation, we developed a stochastic reaction-diffusion system to describe the system in two different scenarios of virus production.In Scenario 1, infected cells produce viruses and DIPs through cell bursting.Therefore the position of virus production is accidental, which leads to a higher probability of patchy formation.In Scenario 2, infected cells keep producing viruses and DIPs.The virus is produced continuously at the cell position, making the spatial distribution concentrated in one point.The patchy pattern observed in the experiments can be regenerated in our stochastic simulation results.Our model provides a good framework for studying reaction-diffusion systems under stochastic effects.
We also built a hybrid algorithm for stochastic simulation.Classical stochastic methodologies are computationally intensive in two-dimensional cases.Our algorithm is based on the pseudo-compartment method [31]  interfaces to adjust complex systems.It combines two scales of simulation methods for modeling the reaction processes and can capture the advantages of the methods with different scales.It improves computational efficiency and maintains critical stochastic features.Our method can provide a numerical framework for studying the spatial stochastic effect of other biological systems and is compatible with different scale stochastic study methods like stochastic differential equations.
We quantitatively studied the spread rate of the virus and showed the relationship between the spread radius growth rate and the initial dose of DIP.To measure the patchiness, we computed the q-statistic.Our simulations can simultaneously capture two spatial spread features (patchiness and spread rate) in wet-lab experimental data, which was not achieved in previous works.It supports that the DIPs slow down the growth of virus particles and make them more patchy.These quantitative methods and statistics are useful tools to understand and explain the diverse spatial-temporal features in complex biological systems.

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

Fig 2 .
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 Ω j P (t) is modeled by PDE and the domain Ω j C (t) is modeled by compartment-based SSA.The amount of each species in a compartment in Ωj P (t) is C k p s (x, t)dx.If any reactant amount is below the threshold θ, then that compartment is part of Ω j C (t).Individuals can move between the boundary compartment of Ω j C (t) and the pseudo-compartment in Ω j P (t).In the two-dimensional case, diffusion takes four directions: up, down, left and right.

3 .
Determine the initial interface for each reaction j, j ∈ {1, 2, • • • , M }: (a) Find all k such that min s∈Sj {N s (k, t)} < θ, where S j contains all species involved in reaction j, then the k-th compartment is part of the stochastic domain Ω j C , and otherwise part of the PDE domain Ω j P .(b) All compartments adjacent to Ω j C (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 −1 , then for the corresponding (s, k), N s (k, t + τ ) = N s (k, t) − 1 and p(x, t + τ ) = p(x, t) + I [x∈C−1] /h.Here, I [x∈A] 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 −1 to stochastic domain and p(x, t) > 1/h for all x 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.

Fig 3 .
Fig 3. A representative example of interpreting the preprocessing of experimental figures.
Fig 4A shows the time series of C * V and C * V D spatial distribution in a 2D domain at time 9-25h with no DIP, as well as images of experimental data at the same time in [3].In both simulations and experiments, viruses are uniformly radially distributed.In Fig 4B, the initial conditions include C V D (0) = 100 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, patchiness is not observed in PDE simulations.March 27, 2023 12/29

Fig 4 .
Fig 4. Dynamics of virus and DIP in cells in PDE simulations and experiments.A: Time series plot of virus in cells (C * V ) and DIP in cells (C * V D ) growth in PDE simulations and representative experiment with initial DIP equal to 0. B: Time series plot of virus in cells (C * V ) and DIP in cells (C * V D ) growth in PDE simulations and representative experiment with initial DIP equal to 84.

Fig 8 shows the radius of virus versus time 9 ≤
t ≤ 25 (h) with different initial DIP inputs in 2 scenarios simulations and experiments.
Fig 9 shows the relationship between the virus radius growth rate and initial DIP dose March 27, 2023

Fig 5 .
Fig 5. Dynamics of virus and DIP in cells in 2 scenarios simulations and representative experiment with initial DIP equal to 0. Row 1 and 2 are time series plots of virus in cells (C * V ) and DIP in cells (C * V D ) growth in Scenario 1 (infected cells produce virus and DIPs through cell bursting); Row 3 and 4 are time series plots of those in Scenario 2 (infected cells keep producing viruses and DIPs); Row 5 is the representative experimental results.

Fig 6 .
Fig 6.Dynamics of virus and DIP in cells in 2 scenarios simulations with C V D22,22 (0) = 4 and representative experiment with initial DIP equal to 1. Row 1, 2 are time series plots of virus in cells (C * V ) and DIP in cells (C * V D ) growth in Scenario 1 (infected cells produce virus and DIPs through cell bursting); Row 3, 4 are time series plots of those in Scenario 2 (infected cells keep producing viruses and DIPs); Row 5 is the representative experimental results.

Fig 7 .
Fig 7. Dynamics of virus and DIP in cells in 2 scenarios simulations with C V D22,22 (0) = 40 and representative experiment with initial DIP equal to 10. Row 1, 2 are time series plots of virus in cells (C * V ) and DIP in cells (C * V D ) growth in Scenario 1 (infected cells produce virus and DIPs through cell bursting); Row 3, 4 are time series plots of those in Scenario 2 (infected cells keep producing viruses and DIPs); Row 5 is the representative experimental results.

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

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

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

Fig A.
Fig A.1.Dynamics of virus and DIP in cells in 2 scenarios simulations with C V D22,22 (0) = 200 and representative experiment with initial DIP equal to 84.Row 1, 2 are time series plots of virus in cells (C * V ) and DIP in cells (C * V D ) growth in Scenario 1 (infected cells produce virus and DIPs through cell bursting); Row 3, 4 are time series plots of those in Scenario 2 (infected cells keep producing viruses and DIPs); Row 5 is the representative experimental results.

1 .
Fig A.1.Dynamics of virus and DIP in cells in 2 scenarios simulations with C V D22,22 (0) = 200 and representative experiment with initial DIP equal to 84.Row 1, 2 are time series plots of virus in cells (C * V ) and DIP in cells (C * V D ) growth in Scenario 1 (infected cells produce virus and DIPs through cell bursting); Row 3, 4 are time series plots of those in Scenario 2 (infected cells keep producing viruses and DIPs); Row 5 is the representative experimental results.

Fig A. 2 .
Fig A.2. Strata used in the computation of q-statistics.

Fig A.
Fig A.3.A representative example to illustrate the image processing of simulation results for computing the q-statistic.A: the spatial distribution of C * V in 2D with the same color bar as Fig 5-7, which is then converted to binary.B: all points where C * V < 50 are black ((0,0,0) in RGB); otherwise, they are red ((255,0,0) in RGB).
Fig A.3.A representative example to illustrate the image processing of simulation results for computing the q-statistic.A: the spatial distribution of C * V in 2D with the same color bar as Fig 5-7, which is then converted to binary.B: all points where C * V < 50 are black ((0,0,0) in RGB); otherwise, they are red ((255,0,0) in RGB).
1. Initialize the time, t = t 0 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.
d) Update the density for the pseudo compartment.Update the density for the pseudo compartment.(c)Updatethecurrenttime, 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 min s∈Sj {N s (k, t)} < θ, where S j contains all species involved in reaction j, then the k-th compartment is part of the stochastic domain Ω j C , and otherwise part of the PDE domain Ω j For the compartment C k that change from PDE domain to stochastic domain, let P s = { C k p s (x, t)dx}.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 (Ω j P ) for x ∈ Ω j P ; • otherwise, N s (k, t I ) = floor (N s (k, t)) and p s (x, t I ) = p s (x, t) + P s / Area (Ω j P ) for x ∈ Ω j P .(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.
(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:

Table 1 .
Parameters used in the simulations.