Mechanistic Modeling of a Novel Oncolytic Virus, V937, to Describe Viral Kinetic and Dynamic Processes Following Intratumoral and Intravenous Administration

V937 is an investigational novel oncolytic non-genetically modified Kuykendall strain of Coxsackievirus A21 which is in clinical development for the treatment of advanced solid tumor malignancies. V937 infects and lyses tumor cells expressing the intercellular adhesion molecule I (ICAM-I) receptor. We integrated in vitro and in vivo data from six different preclinical studies to build a mechanistic model that allowed a quantitative analysis of the biological processes of V937 viral kinetics and dynamics, viral distribution to tumor, and anti-tumor response elicited by V937 in human xenograft models in immunodeficient mice following intratumoral and intravenous administration. Estimates of viral infection and replication which were calculated from in vitro experiments were successfully used to describe the tumor response in vivo under various experimental conditions. Despite the predicted high clearance rate of V937 in systemic circulation (t1/2 = 4.3 min), high viral replication was observed in immunodeficient mice which resulted in tumor shrinkage with both intratumoral and intravenous administration. The described framework represents a step towards the quantitative characterization of viral distribution, replication, and oncolytic effect of a novel oncolytic virus following intratumoral and intravenous administrations in the absence of an immune response. This model may further be expanded to integrate the role of the immune system on viral and tumor dynamics to support the clinical development of oncolytic viruses.


INTRODUCTION
Over the past 4 decades, the oncology treatment landscape has dramatically changed with the development of advanced and targeted immunotherapies, which complement or even replace classical chemotherapy and radiotherapy strategies (Madden, 2018). Oncolytic virus therapy represents one such novel class of immunotherapy. Advantages of this approach rely on the ability of oncolytic viruses to selectively replicate in cancer cells without harming normal tissues. This process leads to a direct lysis of the tumor mass, as well as the generation of a de novo anti-tumor immune response or boosting of an existing one (Kaufman et al., 2015).
At present, talimogene laherparepvec (T-Vec or Imlygic®) with intratumoral (i.t) administration is the only FDAapproved oncolytic viral therapy for the treatment of melanoma patients with injectable but non-resectable lesions in the skin and lymph nodes (Pol et al., 2016). Nonetheless, several oncolytic viruses are currently in clinical development to target solid malignancies (Raja et al., 2018). While i. t. administration maximizes the viral load reaching the tumor, it is only feasible for accessible tumors. On the other hand, the intravenous (i.v.) administration would potentially expand oncolytic viral treatment to less accessible tumors as well as obviate the need for interventional radiology, among other benefits. However, systemic administration needs to overcome viral neutralization, liver and spleen sequestration, and vessel extravasation to exceed the a priori unknown "viremic threshold" above which tumor destruction can be achieved (Russell et al., 2012).
To date, there have been only a few published reports describing i. v. administration of oncolytic viruses in clinical trials. In these studies, viral levels as well as biological activities at tumor lesions were observed in some patients (Kaufman and Bommareddy, 2019), however, no evidence for clinical responses have been observed to date (Breitbach et al., 2011;Russell et al., 2014). Furthermore, the currently available data are from early clinical studies with a limited number of patients wherein the endpoints were focused on safety and tolerability rather than therapeutic activity.
In this regard, prior to embarking on large clinical trials, it may be prudent to understand the biology of oncolytic viruses delivery related to tumor response in more detail (Kaufman and Bommareddy, 2019) and mathematical models can provide a quantitative understanding of the biological processes that an oncolytic virus undergoes following different routes of administration.
Various mathematical models to characterize in vitro viral dynamics (Bagheri et al., 2011;Titze et al., 2017), in vivo behavior limited to tumor size metrics (Bajzer et al., 2008;Eftimie et al., 2011), theoretical frameworks (Wein et al., 2003;Mok et al., 2009;Paiva et al., 2009;Malinzi et al., 2017;Jenner et al., 2018;Cassidy and Humphries, 2019;Al-Tuwairqi et al., 2020) and develop predictive in-silico trials (Jenner et al., 2021) have been so far described. To the best of our knowledge, the integration of kinetic parameters including viral distribution to tumor lesions and comparing i. t. and i. v. administrations has not yet been addressed. A greater understanding of these parameters is critical to optimize dosing, schedules, and routes of administration of oncolytic viruses in clinical trial design.
V937 (formerly named CVA21) is a naturally occurring coxsackievirus currently under clinical development. In vitro (Shafren et al., 1997a) and in vivo (Shafren et al., 2004;Au et al., 2005) oncolytic activity of V937 has been demonstrated in preclinical studies. However, little is known regarding how V937 infection, replication, and tumor distribution in vivo relate to the anti-tumor response. The objective of this work was to develop a mechanistic quantitative framework to elucidate the interplay between viral kinetics, dynamics and viral distribution to tumor lesions and then link this information to an anti-tumor response following i. t. and i. v. administration in human xenograft tumor models in immunodeficient mice. The mechanistic model described here will provide the basis to integrate data related to immune response that would allow elucidation of an oncolytic virus driven anti-tumor response under various treatment regimens including combination with immune checkpoint inhibitors.

Experimental Data
In vitro V937 replication data, in vivo V937 levels in sera and tumors, and tumor sizes in control and V937-treated mice in multiple human melanoma xenograft studies were integrated in the analysis to build the mechanistic model (raw data depicted in Supplementary Figure 1).
Unpublished in-house data and extracted data from published reports were both processed using R (R Project for Statistical Computing, RRID:SCR_001905) v3.6.1 through RStudio interface (RStudio, RRID:SCR_000432) v1.2.1335. In cases where data were only available from graphs, WebPlotDigitizer (WebPlotDigitizer, RRID:SCR_013996) software was used to extract mean profiles over time. Descriptions of the available data are as follows: In Vitro Viral Replication Data of V937 replication in two human melanoma cells lines SK-Mel-28 (NCI-DTP Cat# SK-MEL-28, RRID:CVCL_0526) and ME4405 (RRID:CVCL_C680) and a rhabdomyosarcoma cell line transfected with ICAM-1, RD-ICAM-1, were extracted from Au et al., 2005. In brief, cell monolayers in 6-well plates were treated with approximately 10 6 TCID 50 (half-maximal tissue culture infectious dose) of V937 over 1 h, however the number of seeded cells is not available. At the various time points (0, 2, 4, 6, 8, 10, 12, 24, and 48 h), the cell monolayers were lysed and the viral yield of the cell lysate was determined in triplicate using an end-point titration assay. Average viral titers were digitalized from reported graphs.

In Vivo Viral Kinetics
Viral RNA measurements of V937 in serum and tumor were obtained from an in-house experiment ( Table 1).
Ethical review and approval were not required for the animal study as our work presents a modelling exercise on data already published for which the ethical approval was obtained at the time of experiments were performed (by others in other institutions).
Severe combined immunodeficient (SCID) mice were intradermally inoculated with 2 × 10 6 SK-Mel-28 cells in hind flank on Day −20 for tumor growth. On Day 0, animals were assigned into three groups (n 16/group) with following treatment regimens: 1) a single i. v. dose of PBS (phosphatebuffered saline), 2) a single retro-orbital dose of 10 4 TCID 50 V937 (low dose) 3) a single retro-orbital dose of 10 7 TCID 50 V937 (high dose). Sera and tumors were harvested from two animals per group following euthanasia at various time points. V937 in sera was quantified by both real-time reverse transcriptase polymerase chain reaction (RT-qPCR) (copies/ mL) and cell infectivity (TCID 50 /mL), while V937 in tumor samples was quantified by real-time RT-PCR normalized by total RNA (copies/μg RNA). Viral levels were detected in all samples except for two tumor samples harvested at 3 h (low dose V937), and one sample harvested at 17 days (high dose V937).
In this experiment, tumor size was measured using calipers, and the average tumor volumes per group on Day 0 (n 4), Day 14 (n 4) and Day 17 (n 2) were reported.

In Vivo Tumor Growth Inhibition
In addition to the average tumor volumes obtained from the in vivo V937 viral kinetic experiment, tumor size data from two publications (Shafren et al., 2004;Au et al., 2005) were compiled for the analysis. Experimental designs are detailed in Table 1.
In the studies described by Shafren et al., 2004, data were available from three experimental conditions: 1) single i. t. dose of V937 in mice bearing one tumor lesion, 2) single i. t. dose of V937 in one tumor lesion of mice bearing two tumor lesions, and 3) single i. v. dose of V937 in mice bearing two tumor lesions. In all cases, non-obese diabetic (NOD) SCID female mice were subcutaneously (s.c.) injected with 2 × 10 5 SK-Mel-28 cells in one or two flanks, and single injection of phosphate buffered saline (PBS) or V937 (10-10 5 TCID 50 ) was administered when tumors reached a predefined volume. Tumor volume was calculated using the formula for a spheroid. Tumor volumes, were recorded at regular intervals, and average tumor volume were digitalized for analysis.
In vivo growth rate of SK-Mel-28 in NOD-SCID mice following implant of 10 6 cells in one flank from Au et al., 2005 was combined with control groups in SK-Mel-28 study from Shafren et al., 2004 and used for model building.
In addition to SK-Mel-28 cell line, data from control and treated animals inoculated with ME4405 cells were available from Au et al., 2005. This information was used as validation (see Experimental data section of Supplementary Methods and Results for detailed description of experimental design).

Mechanistic-Based Model Building
The mechanistic model was built following a sequential and integrative approach (Figure 1). First, in vitro data were used to characterize dynamics of viral infectivity and replication. Following, in vivo viral kinetics, viral distribution in tumor and tumor growth were described taking into account the viral dynamic model previously developed.
The model, depicted in Figure 2, accounted for the following biological processes: 1) viral clearance from sera 2) viral distribution to tumor mass, 3) proliferation of tumor cells 4) uninfected tumor cells (uTC) to be infected by the virus, 5) viral replication in infected tumor cells (iTC), and 6) viral induced death of infected tumor cells.
The following subsections describe, 1) the final model structure and parameter estimates, 2) data analyses as methodology used to build the model, 3) model selection and evaluation.
Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 705443 2002b), the model described the time course of three main entities: uTC, iTC and viral load in the extracellular medium (VL EC ) as shown in eqs. 1-3. uTC can proliferate with a net growth first order rate constant (λ) and be transformed to iTC in the presence of virus in the medium. This infectivity process takes place through a second order rate constant (β) that consumes both, viruses and uninfected tumor cells. Once cells are infected, V937 can replicate within iTC and generate copies of virions (α) per infected cells during its life span. Subsequently, iTC will be killed by viral induced cell death at a rate (δ). Extracellular virions can also be degraded with a first order rate constant (K VIR ).
The system was initialized to the administered dose (TCID 50 ) in the compartment of VL EC , and number of seeded cells at confluence for uTC (1.2 × 10 6 cells/well 1 ). iTC were considered to be zero at baseline.
Total predicted viral load, which is computed as the sum of viral load in the medium (VL EC ) and inside cells (α × iTS), was fitted to measured viral concentrations from the in vitro replication experiment, using the volume of 3 ml per well of a 6-well plate.
Given the short duration of the in vitro evaluation, viral degradation and cell proliferation processes were considered FIGURE 1 | Schematic representation of the modelling and data workflow, highlighting the key processes identified at each step.

In Vivo Viral Kinetic, Viral Dynamics and Tumor Growth Inhibition Model
To characterize V937 viral kinetics in vivo, the time courses of V937 viral load in sera (VL S ) and tumor vasculature (VL VAS ) were modeled through a minimum physiological model (eqs. 4-5). Mouse-specific tumor blood flow (Q TUM ;144 ml/d) (Baxter et al., 1995) and real organ volumes were used in the analysis. Tumor vascular volume (V VAS ) was assumed to be a 7% of total tumor volume at baseline and serum volume (V S 0.974 ml) was fixed at 55% (Baxter et al., 1995) of blood volume in mice [77 ml/kg for a 23 g mouse (Mitruka and Rawnsley, 1981)]. Viral elimination rate (K VIR ) and tumor retention factor representing viral affinity for the tumor (RF) were estimated as drug-specific parameters.
The initial values of VL S and VL VAS at time 0 were set to the administered dose (TCID 50 ) and 0, respectively, following i. v. dose. In the case of i. t. administration, VL S was initialized to 0 and VL VAS to the administered dose at time 0. The kinetic model was linked to the previously developed viral dynamic model at tumor level replacing VL EC by VL VAS in eqs. 1-2 as follows: As for the in vitro model, uTC compartment was initialized to the number of initial tumor cells that were derived from measured tumor size at baseline, and iTC was assumed to be zero at baseline.
To identify K VIR and RF drug parameters, predicted viral load concentration in serum (VL S /V S in TCID 50 /mL or copies/mL after correction by a factor accounting for the ratio between TCI50 and copies, referred to as RATIO) and in tumor [computed as the sum of vascular, VL VAS , and intracellular levels, α × iTS, normalized by an estimated factor to account for the amount of RNA (µg) per tumor volume unit, RNAVOL], were fitted to experimental levels of V937 in sera and tumor over time (see In Vivo Viral Kinetic Experimental Data section). The parameters α and β were fixed to those obtained in vitro, while λ and δ parameters, considered to vary between in vitro and in vivo experimental scenarios, were estimated by fitting tumor size model predictions (uTC + iTC) to experimental tumor size data (see Tumor Growth Inhibition Experimental Data section). Note that tumor size was measured in mm 3 while uTC and iTC are expressed as cell units, therefore a standard conversion factor of 10 6 cells per mm 3 was used (Makkat et al., 2007). In addition, a proportionality ratio (RATIO) between TCID 50 and copies was computed to account for the conversion between both metrics as the model describes the dynamics of TCID50/mL and copies/mL, assuming that both variables have same time course but different magnitude.

Data Analysis
For in vitro or in vivo viral levels, each observation was obtained from a different well or animal; therefore, a naïve pool approach was used (i.e. data were pooled together and analyzed as coming from one single ID). For tumor size data in vivo, each mean profile was considered as an independent animal, thus enabling the use of the population approach to characterize inter-subject variability in model parameters as well as residual error.
The software NONMEM 7.4 (NONMEM, RRID: SCR_016986) and the First Order Conditional Estimation method with interaction algorithm were used for the analysis.
Data were logarithmically transformed during the analysis and an additive error model in the logarithmic scale was used to account for the discrepancies between model predictions and observations. A different residual error for each measurement was considered. When the population approach was applied the inter-subject variability was modelled exponentially (Kiang et al., 2012;Mould and Upton, 2012).

Model Selection and Evaluation
Model selection was largely driven by biological plausibility and capability of the model to describe the tendency of the data, taking into account also 1) parameter precision, 2) classical goodness-of-fit plots including observation versus model prediction and conditional weighted residuals over time or over model predictions, and 3) objective function value (OFV, approximately equal to -2xlog-likelihood). When models were nested, a drop of 3.84 or 6.63 in OFV was considered significant at the levels of 5 and 1%, respectively.

Model Exploration
The dynamics of the different entities of the system were explored over the range of evaluated dose levels (10-10 7 TCID 50 ), comparing outcome after i. v. or i. t. administration including one or two tumor lesions, the latter allowing for exploration of abscopal effects.
In addition, a parameter scan was performed varying one parameter at a time to explore its impact on viral levels and tumor response for selected design scenarios. During the simulation step,, a function (FMOI) was introduced at the infectivity term: FMOI × β × uTC × VL VAS ; this was used to avoid infectivity of cells at very low virus levels just because high number of uninfected tumor cells in a continuous manner, rather than setting viral levels to zero at random value. This function cancels the infection of new cells when the viruses are too low compared to the number of total cells, with FMOI defined as follows: Where MOI (multiplicity of infection) is computed as the ratio between viral load levels in the vasculature and total tumor burden (iTS + uTS) and MOI 50 represents the MOI at which 50% of maximum infectivity is obtained. After a sensitivity analysis, MOI 50 was fixed to a low value (10 -6 ) that did not affect model characterization of the experimental data.

Sensitivity Analysis
To evaluate the impact of simultaneously varying all parameters on model output, a global sensitivity analysis was performed following the approximation technique described by Saltelli et al., 2008;Saltelli et al., bib_Saltelli_et_al_20102010 of the Sobol's method (Sobol′, 2001) using the SAFE toolbox available as R package (Pianosi et al., 2015), where to have even probabilities of sampling parameters differing in several order of magnitudes the latin hypercubic sampling method from a uniform distribution was used. Tumor size at day 14 was selected as the variable to represent drug response for model output.

Model Applicability
The combined impact of the most influential parameters was explored for two-by-two combinations of different dose levels and dosing routes. To do so, a virtual population (n 1,000) was simulated varying all other parameters but those of interest (i.e. α and β, α and RF or β vs RF) following the same sampling approach as described for the global sensitivity analysis.
Tumor size profiles of the virtual population were then simulated using combinations of the two parameters of interest over the plausible parameter space. Probability of response under the different scenarios (i.e. combination of parameters, dose levels and dosing routes) was computed as the probability of observing at least 20% tumor shrinkage at day 14.

In Vitro Viral Dynamic Model
The levels of viral load obtained in vitro showed an initial decay reflecting infectivity, followed by a rapid increase up to 10 8 TCID 50 /mL, 8-12 h post-infection, indicating viral production (raw data depicted in Supplementary Figure 1A).
Due to the availability of only total viral levels, it was not feasible to account for intracellular viral production and subsequent release to the cytoplasm as independent processes. Describing viral replication as a unique process dependent on both cell death and number of copies per cell allowed us to 1) obtain precise estimates of viral infectivity and production and 2) account for intracellular levels without increasing model complexity. As expected from the quick viral plateau reached in vitro, a fast death rate constant (δ) was identified, although with low precision due to the lack of cell death time profiles. Accounting for degradation of virus in the medium provided a very low estimate of viral clearance (K VIR ), suggesting this process could be neglected without affecting model performance Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 705443 (p > 0.05), probably as a consequence of the short duration of the experiment. Different α, β or δ parameter estimates for the three available cell lines (RD-ICAM-1, SK-Mel-28 and ME4405) were evaluated. Best results in terms of model performance and parameter precision were obtained using a different δ, estimated to be higher for ME4405 than for RD-ICAM-1 and SK-Mel-28. The use of the different δs provided a significant drop in objective function (p < 0.05) and improved parameter precision (from 100% error to less than 60%) ( Table 2). Nonetheless, α and β estimates did not differ statistically across the three cell lines. Overall, a good description of the in vitro data was obtained ( Figure 3A).

In Vivo Viral Kinetic and Tumor Growth Inhibition Model
I.V. administration of V937 at both 10 4 and 10 7 TCID 50 resulted in a rapid increase in V937 in sera and tumors, 6-24 h post administration. Despite the 3 log10 difference in dose (10 4 vs 10 7 TCID 50 ), similar viral concentrations were observed in plateau suggesting viral capability to replicate in vivo (raw data depicted in Supplementary Figure 1B).
A semi-physiological pharmacokinetic model that accounts for viral replication in the tumors provided a sufficient description of the available kinetic data in both serum and tumor. Infectivity (β) and replication (α) estimates were fixed to those obtained from the in vitro analysis, while physiological parameters (i.e. tumor blood flow, serum volume and fraction of tumor vasculature) were taken from the published literature. The rest of the model parameters were estimated. High V937 elimination (t 1/2_KIV ∼4.3 min) and V937 tumor retention (RF), supported by the low initial viral levels and the capability of the virus to efficiently distribute to tumor and quickly replicate, were identified with good precision.
Regarding tumor size data, certain dose response could be identified, but no major differences across routes were experimentally detected, though large variability in tumor size at dosing time was observed (raw data depicted in Supplementary Figure 1C). When simultaneously modelling control and treated groups, exponential growth provided a better overall description (p < 0.001). Variability across mean tumor profiles was identified on baseline tumor cells and death rate of infected tumor cells. In addition, estimating different baseline and growth rate for secondary tumor lesions produced a significant (p < 0.01) decrease in OFV.
Acknowledging that tumor growth/shrinkage could also impact V937 viral distribution to the tumor (i.e. tumor blood flow), different models that assumed changes in Q TUM with tumor mass were explored either assuming a proportional relationship between them (i.e. tumor blood flow and tumor mass) or as a power model that allow to account for the fact that not all tumor mass is perfusable (Ferl et al., 2005). None of the models was supported by the data, neither provided a different description of pharmacokinetic (PK) profiles and therefore this feature in the model was not included.
Overall, the final model was able to satisfactorily describe all sources of in vivo data simultaneously ( Figures 3A-C) enabling precise parameter estimates (Table 2). Moreover, the final model structure could adequately describe the limited experimental tumor response data obtained when implanting ME4405 cell line in xenografts (see results section of supplementary material).

Model Exploration
The final model was used to explore the behavior of the observed and unobserved elements of the system under different experimental scenarios. Following i. v. or i. t. administration, V937 reaches a quick equilibrium of distribution between serum (VL S ) and tumor vasculature (VL VAS ) as shown in Figure 4 and at both at the primary and secondary lesion as shown in Supplementary  Figure 2. Given the estimated high tumor affinity and infectivity, the model predicts that the vast majority of cells get quickly infected (in less than 1 h at a dose of 10 4 TCID 50 ) and start to produce new viral copies that return to serum. This aspect, observed at all experimental doses available ( Supplementary  Figure 3), explains the increment in observed circulating levels, despite the fast viral degradation. Time to reach the maximum infectivity of tumor cells, and thus peak viral load levels, depends on both dose and initial tumor burden, with a shortening in time observed as dose level (Supplementary Figure 3) or tumor size increases (data not shown). Under these conditions, differences between routes of administration are mainly located at early time points, before maximum infectivity is achieved, but in all cases tumor response is ultimately attained, as also observed experimentally.

Sensitivity Analysis
The impact of varying one model parameter at a time within a plausible parameter space was evaluated as exemplified for α in Figure 5A and for the rest of model parameters in Supplementary Figure 4. Only changes in viral replication (α), viral infectivity (β) or tumor retention factor (RF) were able to invert the course of tumor response from cure to progression when parameter values were decreased. The rest of parameters showed an impact on the rate of response, like the infected cell death rate (δ), or minor influence on overall profiles, without changing tumor progression.
When performing a global sensitivity analysis to simultaneously explore the impact of model parameter changes on tumor response, similar results were obtained ( Figure 5B). α, β, RF and δ were the most influential parameters, accounting for around 10-16% of the variability in response as single parameters (first-order index), and 28-40% in total when taking into account the different doses evaluated and routes of administration (Supplementary Figure 4). These data highlight the importance of combined effects among parameters to explain model outcome. However, a different scenario is observed for the lowest dose of 10 TCID 50 ; this is due to the lack of infectivity and response for most of the simulated parameter vectors, thus leaving the rate of tumor growth and baseline tumor size as the most influential parameters determining tumor size outcome.

Model Applicability
The model was further evaluated to study the probability of observing a response (i.e. tumor shrinkage greater than 20%) for different two-by-two combinations of most influential model parameters with an impact on ultimate response. This methodology allowed us to identify those areas where the probability of observing a response is almost null or very high, regardless of the combination of all other model parameters. As expected under this system, similar profiles were obtained after i. v. and i. t. administration, with slightly higher probability of response for the latter (data not shown). In all explored scenarios, increasing the dose translated into an increase in the probability of response, covering a wider area of the parameter space ( Figure 6).

DISCUSSION
Mathematical modelling in immunology and virology has been extensively explored to provide a better understanding of the infection time course and to optimize viral therapies (Perelson, 2002a;Santiago et al., 2017). With the increasing interest in the use of oncolytic viruses in treating cancer, different models that integrate these concepts into a tumor dynamic environment have appeared (Santiago et al., 2017).
In this work, we present a modelling framework using ordinary differential equations that is capable of simultaneously accounting for three key processes of a novel oncolytic virus, V937, currently in clinical development: 1) capability to infect and replicate within tumor cells, 2) pharmacokinetics including distribution of V937 virus to tumor and 3) finally oncolytic V937 effects in vivo. To the best of our knowledge, this is the first time that distribution to tumor has been explicitly taken into account to enable characterization of different routes of administration.
To build the framework, an intermediate approach between the classical top-down and bottom-up approaches (i.e. data driven vs biology driven) was followed, thus balancing mechanistic knowledge and available data. Data from different sources, levels of information and even measurement units were integrated, highlighting the importance of quantitative models to combine and exploit experimental data as already shown in other therapeutic scenarios (Campagne et al., 2018;Parra-Guillen et al., 2020).
In this respect, the use of in vitro data to identify viral infectivity and replication processes, which are hard to identify only using in vivo data, has been of particular relevance. Parameter comparison with other viruses is not straight-forward as different dosing and measurement units (pfu, copies, TCID 50 ) are used across experiments. However, in absolute terms, V937 replication capability was within the reported ranges for oncolytic viruses of 50-1,350 virions/cell (Liu et al., 2000), which in most cases come from direct experimental observations and not from modelling exercises. Similarly, estimated in vitro death rate was in line with the upper ranges reported by Titze et al., 2017 (1.2-761 1/h) using a similar model structure, as well as the infectivity rate (0.35 × 10 −8 to 0.98 × 10 −8 ). Nonetheless, larger intervals depending on the virus and cell type can be found in the literature for viral infectivity (10 −7 -10 -10 ) (Bajzer et al., 2008;Mahasa et al., 2017;Cao et al., 2018).
V937 employs the intercellular adhesion molecule I (ICAM-I, CD54) receptor for attachment and viral entry (Shafren et al., 1997a;1997b); this receptor is overexpressed in numerous malignant cells, including melanoma (Kageshita et al., 1993;Hayes and Seigel, 2009). V937 can also bind to the DAF receptor; however, it appears that DAF may function as a lowaffinity attachment receptor either enhancing viral presentation or providing a viral sequestration site for subsequent high-affinity binding to ICAM-1 (Shafren et al., 1997b). Despite different ICAM-I expression in the three in vitro cell lines, differences in viral entry across cell lines could not be adequately identified. This could potentially be due to the experimental data and design, e.g. the lack of dose range as well as errors associated with end point titration for TCID 50 assessment. However, it could also be due to a sufficiently high ICAM-1 expression levels in all cases (Au et al., 2005) compared to the threshold of 5,000 ICAM-1 molecules per cell required for V937 to exhibit in vitro activity (Annels et al., 2018). Nonetheless, and given the mechanistic nature of the model, additional data for example regarding ICAM-I expression could be easily introduced to further refine the model.  In the current model tumor cells infection is predicted to take place quickly and in less than a day due to either the high infectivity rate constant obtained in vitro in combination with the high predicted tumor levels after i. t. administration, or the fast tumor distribution and high tumor retention estimated after i. v. distribution. That critical aspect of the model has been explored carefully and is supported by the data. Specifically, viral replication in vitro indicated a fast infectivity process, with increasing viral levels reaching a plateau 8-12 h after V937 exposure. Moreover, it has been described that the virus is capable to trigger an extensive lytic destruction of melanoma cells 23 h after infection at a dose level of one TCID50/cell (Shafren et al., 2004), level similar to the one used in vitro. Regarding the vivo scenario, a sharp and sustained increase in viral levels at tumor is also observed after i. v. administration, despite the high viral clearance, supporting that infection and replication of the virus takes place fast, given the current model structure.
One novel aspect introduced in this work is the description of viral kinetics and viral distribution to tumor using a semiphysiological pharmacokinetic model that enables simultaneous description of both i. v. and i. t. administration routes. Viral measurements at tumor level are scarce in the literature (Kelly et al., 2008;Workenhe et al., 2014;Garcia-Carbonero et al., 2017;Andtbacka et al., 2019), and even more so over time. However, these types of data have been essential to characterize the tumor compartment in the model. In our modelling framework, the tumor was represented by two subcompartments: 1) a vascular compartment in fast equilibrium with serum and 2) a cellular compartment which gets infected upon viruses arriving to the vascular compartment. A high retention of the virus at tumor level (RF 623), which cannot be compared due to lack of literature data, and direct release of newly formed virions to the vascular compartment were required to successfully characterize the data behavior, suggesting no distribution limitations between serum and tumor in our preclinical setting. Moreover, an estimated half-life of 4.3 min, in line with the fast clearance observed for oncolytic viruses in clinical studies (Garcia-Carbonero et al., 2017), was obtained. This result indicates that the sustained levels observed in tumor and plasma are due to V937 capability to infect and replicate in immunodeficient mice, which ultimately leads to a potent oncolytic response as reflected by the infected cell death rate constant (δ) estimate, which is in line with those reported in the literature (Okamoto et al., 2014;Cao et al., 2018).
One of the major hurdles in current drug development in oncolytic viruses, and immune-therapies in general, is the need for predictive animal models (Russell et al., 2012). Certainly, xenograft mouse models that lack an immune system can be limited. However, they present a valuable opportunity to assess the properties of the oncolytic viruses in isolation, without any other limitation, on human tumor cell lines, as illustrated in this work. In this mouse model, viral retention, infectivity and replication at tumor level were identified as the key processes controlling V937 tumor response. These processes are tightly interconnected and difficult to identify simply using tumor response data. Viral infectivity and replication identified in vitro experiments could be directly integrated into the in vivo model structure to enable an adequate prediction of tumor response, thus providing an in vitro/in vivo framework for oncolytic viruses that can be used to support the selection between candidates based on their in vitro properties.
In this tumor mouse model and despite the rapid systemic clearance, minor differences between i. v. and i. t. routes of administrations were observed due likely to high viral infectivity and replication in the tumor level, as reflected by the corresponding parameter estimates, which in principle can lead to complete tumor eradication. However, this result should be interpreted with caution at the time to translate to other scenarios as the impact of the immune system is not yet considered and its role on the viral infection and proliferation mechanisms cannot be ruled out. In addition, given the small number of animals used in the experiments further data would be needed to validate the current model structure and parameter estimates. In this regard, model development should be seen as an iterative process that needs to be coupled with experimental work in order to reflect biology of the system and maximize model usefulness. As a next step, information from syngeneic mouse models that include not only tumor size measurements, but also kinetic levels in serum (and ideally in tumor as well) and relevant immune response markers, could be integrated into this framework. Such developments could facilitate a mechanistic and quantitative understanding of the dual role that the immune system can play in viral response, potentially limiting viral infectivity as well as triggering a potent anti-tumor immune response.
In summary, a mechanistic framework integrating in vitro viral dynamic properties into an in vivo system to describe oncolytic effects of V937 on tumor response in immunodeficient mice has been successfully developed. This model allows for a better understanding of the role that the different processes play in the final outcome, and enables selection between oncolytic virus candidates based on in vitro/ in vivo features, such as infectivity. Moreover, the developed model can serve as a backbone to include future additional biological components, such as the immune response, to provide a quantitative understanding of the balance between immune and antiviral response. This would facilitate a better understanding of the limitations of systemic administration in immunocompetent scenarios, guide dosing strategies, and help identify potential combination strategies to ultimately support the development of programs for oncolytic viruses.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
Ethical review and approval were not required for the animal study as our work presents a modelling exercise on data already