# Computron: Serving Distributed Deep Learning Models
with Model Parallel Swapping

Daniel Zou${}^{\text{1, *}}$ Xinchen Jin${}^{\text{2}}$ Xueyang
Yu${}^{\text{2}}$ Hao Zhang${}^{\text{3}}$ James Demmel${}^{\text{1}}$
${}^{\text{1}}$UC Berkeley ${}^{\text{2}}$ShanghaiTech University
${}^{\text{3}}$UC San Diego

###### Abstract

Many of the most performant deep learning models today in fields like language
and image understanding are fine-tuned models that contain billions of
parameters. In anticipation of workloads that involve serving many of such
large models to handle different tasks, we develop Computron, a system that
uses memory swapping to serve multiple distributed models on a shared GPU
cluster. Computron implements a model parallel swapping design that takes
advantage of the aggregate CPU-GPU link bandwidth of a cluster to speed up
model parameter transfers. This design makes swapping large models feasible
and can improve resource utilization. We demonstrate that Computron
successfully parallelizes model swapping on multiple GPUs, and we test it on
randomized workloads to show how it can tolerate real world variability
factors like burstiness and skewed request rates. Computron’s source code is
available at https://github.com/dlzou/computron.

<EMAIL_ADDRESS>
##  1 Introduction

In recent years, researchers and practitioners have dramatically improved the
performance of deep learning models, particularly large language models
(LLMs), using two techniques: massive parameterization and fine-tuning. Many
pre-trained models with billions of parameters have been released, and each of
them is being customized for a myriad of tasks through fine-tuned variants. In
a plausible scenario, organizations would host many of these large models,
each similar in architecture and size but tuned to some specific task, to
serve the needs of their internal personnel and external users.

The usual way to serve a model using GPUs is to keep all of its parameters in
GPU memory so that inference runs directly on the accelerator device. When a
model is too large to fit in a single GPU’s memory, a common technique is to
distribute it to multiple GPUs through model parallelism. The amount of memory
onboard each GPU is limited, so an organization would need to purchase many
GPUs to serve all of its models, which could be quite expensive. Worse, the
costly hardware is underutilized when some models receive requests at low or
irregular rates.

Among existing ML serving systems, some such as AlpaServe [7] and Energon-AI
[4] employ model parallelism to serve large models in a distributed fashion.
There are also systems like Clockwork [5] that use memory swapping to overcome
the limitation of GPU memory and improve utilization.

In this paper, we present Computron, a prototype serving system that unifies
model parallelism and swapping. Computron makes it possible to serve multiple
large distributed models that, in total, can exceed the memory capacity of the
GPU cluster they share. In terms of usability, Computron supports Colossal-
AI’s [6] functionality for easy model parallelism during model development,
and it integrates with asynchronous Python frameworks for service deployment.
We discuss several ordering and synchronization problems that constrain the
design of such a system, and we explain how our design solves these problems
to achieve parallelized swapping on multiple GPUs. We evaluate Computron in
two ways. First, we isolate the swapping component to demonstrate that model
parallel swapping does in fact reduce the time taken to load a distributed
model into GPU memory. Second, we test Computron under more realistic
conditions on randomly generated workloads that simulate conditions where
requests may be bursty and skewed to a subset of models.

##  2 Background and Problem

Deep learning is being rapidly adopted in countless business and scientific
applications. In many of these applications where some form of service
deployment is involved, the serving system is a crucial component of the deep
learning workflow. These systems generally operate in a request-response
paradigm by listening to inference requests, running the requested model on
inputs on specialized hardware such as GPUs, then responding with the output.
The design of such a system revolves around the tradeoff of reducing latency
so that end users experience less waiting versus increasing efficiency to save
operational cost, all without compromising model accuracy.

Recent research in language models has popularized the practice of fine-
tuning, where either a portion or all of a pre-trained model’s parameters are
trained on new data from a specific task in order to achieve higher accuracy.
For example, the pre-trained BERT [3] model can be fine-tuned to a variety of
language understanding tasks—from text classification to part-of-speech
tagging to natural language inference—just by retraining the last layers, and
GPT-3 [2] has been fine-tuned on human feedback so that the resulting
InstructGPT [10] model can better align to user intentions. Fine-tuning
generally involves none or minor modifications of the model architecture. As
fine-tuned models become commonplace, so will workloads that involve serving
multiple models with highly similar memory footprints and access patterns.

A second trend spearheaded by the language modeling community has been to
increase model size in pursuit of better accuracy and generalizability. At the
extreme, Megatron-Turing NLG [11] contains 530 billion parameters, and the
comparatively space-efficient LLaMA [12] still contains up to 65 billion
parameters. On Chatbot Arena [14], a platform where humans rate the quality of
language model outputs, most models with competitive performance have at least
6 billion parameters.

Serving multiple instances of such large models would exceed the memory
capacity of all consumer GPUs and many high-end ones. The standard solution is
to distribute large models across multiple GPUs using model parallelism. Two
forms, tensor parallelism (TP) and pipeline parallelism (PP), are well-studied
and commonly used in training workloads [8], but they are just recently
beginning to see use in serving systems as well. Even when all models can fit
into aggregate GPU memory, the AlpaServe [7] team found that there is still
reason to use model parallelism in a serving system because it can reduce
latency in real world workloads.

Another challenge of serving deep learning models is that real world serving
workloads are often unpredictable. Request arrival distributions may be
bursty. Furthermore, across multiple models, request rates may be skewed—some
models may receive a lot more requests than others—and the rates may shift
over time. Hosting all of these models in GPU memory leads to underutilization
in the face of dynamic request patterns, as resources are over-provisioned to
models with low request rates. The resource allocation imbalance is worse for
localized services that expect irregular traffic, and the hardware cost is
higher when larger models distributed across multiple devices are involved.
Both of these factors act against the interests of smaller organizations.

We survey a number of prior works and find that while there are many designs
from which to take inspiration, to the best of our knowledge, no single system
addresses both the model parallelism and the resource utilization challenges
that we have outlined. AlpaServe [7] and Energon-AI [4] are capable of serving
large models using model parallelism, but they host all models on GPUs
following some form of static assignment and are thus bounded by the amount of
available GPU memory. Clockwork [5] serves many deep learning models on a
limited number of GPUs by swapping models between CPU and GPU memory, which
works well when the models are on the order of millions of parameters or less.
However, this approach is not suited for larger models with billions of
parameters that can take several seconds to transfer. ZeRO-Inference [1]
parallelizes parameter transfers across GPUs in order to multiply CPU-GPU
bandwidth, but it is meant for individual layers within a single massive
model.

##  3 Design

In §2, we provide motivation for the problem of serving multiple distributed
deep learning models, and we identify the key issue of resource
underutilization when handling bursty, skewed requests. To deal with these
challenges, we propose a serving system that co-locates multiple models on the
same cluster of GPUs and dynamically swaps distributed model parameters
between CPU and GPU memories. Active models are swapped into GPU memory so
that requests can be served quickly, while unused models are swapped out to
reduce the unnecessary consumption of resources.

The underlying technique of offloading unused models to CPU memory has already
been applied by state-of-the-art serving systems like Clockwork [5] to great
effect. Our particular approach is comparable to demand paging in the context
of virtual memory management; when a model is requested whose parameters do
not currently reside in GPU memory, a replacement policy is used to pick
another model to swap out, and then the requested model is swapped in. Like
other systems, we assume the existence of large CPU memory to hold unused
model parameters. More sophisticated fetching algorithms can be used here
instead, but are beyond the scope of this paper.

Inspired by prior works, we seek to investigate whether model parallelism can
also be used to reduce the latency of model swapping. We hypothesize that on
systems where GPUs have independent PCIe links to the CPU, by specifying
higher degrees of TP and PP to distribute model shards across more GPUs, model
parameter shards can be loaded in parallel to take advantage of greater
aggregate link bandwidth between the CPU and GPUs. A similar optimization is
used by ZeRO-Inference. Should this prove true in practice, it would become
feasible to perform dynamic swapping while serving large models that share a
group of devices, just like smaller models. However, a number of design
considerations arise when we put our hypothesis to the test.

###  3.1 Architecture

Figure 1: Computron architecture for $TP=2$, $PP=2$. One worker is launched
per GPU. Two models, A and B, are co-located in the same parallel
configuration.

Our system uses an engine-worker architecture to manage multiple distributed
model instances at the same time. The centralized engine receives, queues, and
waits for the completion of all model requests. Workers are launched per GPU
in accordance with a user-provided parallel configuration (TP and PP
dimensions) to manage shards of model parameters. Because we assume all models
have a similar size and architecture, we make the simplification to co-locate
each distributed model instance onto the GPU cluster using the same
configuration. Fig. 1 gives an example of the architecture for a $TP=2$,
$PP=2$ configuration.

When the engine receives a request for some model, it pushes the request
object along with a timestamp into a queue specifically for that model.
Concurrently, the engine repeatedly picks a queue to pop oldest request
objects, then packs and submits them to workers in the first pipeline stage as
a single batch entry. Workers at each pipeline stage evaluate batch entries in
submitted order, up to the last stage, at which point last stage workers send
the batch output back to the engine. PP communication occurs through FIFO
pipes, while TP communication is done through distributed collectives, as
represented by the arrows in Fig. 1.

###  3.2 Model Parallel Swapping

Figure 2: Broadcasted load entry violates load dependency. Figure 3:
Synchronous load entry reduces loading parallelism and causes unnecessary
blocking.

In our design, the responsibility of making swapping decisions is delegated to
the engine, so in addition to submitting batch entries, the engine can
initiate another type of action through what we refer to as load entries. A
load entry commands a worker to either load or offload the parameters of an
instance.

Challenges arise when designing how load entries should be submitted to and
processed by distributed workers. A model can only be evaluated on a batch
entry after the model’s parameters are loaded into GPU memory, so as the
engine schedules batch and load entries in some order it deems correct,
workers must respect the load dependencies of that schedule. Furthermore, data
dependencies between adjacent pipeline stage workers delay when later stage
workers receive batch entries, ruling out certain designs like broadcasting
the load entry, as illustrated in Fig. 2. These load and data dependencies are
resolved if workers synchronously process load entries in pipeline order just
like batch entries, but this naive solution has two issues: a batch entry to
some model is unnecessarily blocked by load entries to another unrelated
model, and no loading parallelism is achieved by workers of different stages
in the same pipeline, as shown in Fig. 3.

Figure 4: Comparison of how batch entries and load entries are processed in a
linear worker pipeline.

We propose an asynchronous mechanism for handling load entries that mitigates
these issues. After being submitted by the engine, load entries are pipelined
through worker stages just like batch entries, but a worker does not wait for
loading to complete before passing the load entry forward to the next stage.
This can be done using the stream feature of the CUDA programming model. On
top of the default CUDA stream that executes kernels for model inference, each
worker launches two additional streams to run loading and offloading
operations concurrently. A load entry is completed when every worker finishes
loading/offloading and sends a response back to the engine. The engine is
responsible for avoiding load dependency violations by making sure batch
entries for a model are submitted to workers only after that model has been
fully loaded. This design allows a later batch entry to proceed without
waiting for a previous load entry involving another model to complete, and
this also enables workers of different stages to load shards of a model’s
parameters in parallel. The paths of batch and load entries through one branch
of the system pipeline are depicted in Fig. 4.

One more detail is the use of pinned memory. CUDA requires the CPU-side data
buffer to be in page-locked memory during CPU-GPU data transfers to prevent
interruptions caused by paging. Data objects on CPU are stored in paged memory
by default, so data transfers would incur an extra copy on the CPU side from
paged memory to page-locked memory. We eliminate this extra data movement by
making sure that when a model is offloaded, the parameters are kept pinned in
CPU memory.

##  4 Implementation

We build Computron as a serving system that supports model parallel swapping
based on the considerations presented in §3. We borrow some components from
Energon-AI [4], such as the RPC-based FIFO pipe implementation used for
communication between pipelined worker stages. Just like Energon-AI, Computron
is compatible with Colossal-AI [6] functionality, meaning that users can
easily incorporate model parallelism in their models with minimal changes to
PyTorch source code.

As Computron launches, Colossal-AI automatically handles setting up the
context and communication groups for model parallelism, and it does so using
the same configuration for each instance. The engine is implemented using
Python’s asyncio library, and request scheduling is done in a completely
asynchronous fashion. Because of this, Computron integrates with asynchronous
Python web frameworks such as FastAPI. Requests are scheduled in batches based
on the oldest timestamp, and model swapping uses an LRU replacement policy.

##  5 Evaluation

We design two separate sets of experiments in order to characterize
Computron’s performance. In the first set of experiments, we intentionally
induce the worst case for handling each request and measure how the time to
swap models scales with model parallelism. In the second set, we generate
simulated request workloads using a random arrival process to study how
Computron handles more realistic scenarios.

We conduct experiments on a single GPU node of the Perlmutter supercomputer
managed by NERSC. The GPU node has one AMD EPYC 7763 CPU and four NVIDIA A100
GPUs, each connected to the CPU through a PCIe 4.0 x16 link [9].

###  5.1 Swapping Latency

In §3, we hypothesize that model parallelism linearly decreases the time taken
to load and offload model parameters between CPU and GPU memory. To check this
hypothesis, we design an experiment that forces the worst case scenario where
each request must perform a swap. We launch models concurrently and configure
the engine to only allow one model to reside in GPU memory at any given time.
We then send alternating blocking requests to the two models and measure the
times taken to swap and execute a model at each request. The model size is
fixed in order to test the strong scaling properties of CPU-GPU swapping. The
model chosen is OPT-13B [13], an open-source pre-trained transformer language
model released by Meta AI. Using half-precision floats, OPT-13B has a memory
footprint of about 24 GB. This model is chosen because it can fit into the
memory of a single A100 GPU, which serves as a baseline for swapping time.

Before running the experiment, we estimate the lower bound for how long
swapping should take for comparison. Each CPU-GPU link has a bandwidth of 32
GB/s, so a single GPU is expected to load or offload an OPT-13B model instance
in $24/32=0.75$ seconds. On our test system, aggregate CPU-GPU bandwidth
increases linearly with the number of GPUs, so as the model is distributed to
more GPUs using either TP or PP, the swapping time is expected to inversely
decrease. Swapping time includes both the offloading of one model and the
loading of another, and because our asynchronous implementation overlaps the
two, we measure from when the offload entry is submitted to when both offload
and load entries are completed.

Figure 5: Swapping latency with changing TP scale.

With the theoretical lower bound in mind, we first run the experiment with
three trials that scale the degree of model TP. We use a small input token
length of 2. Fig. 5 visualizes the results of these trials on $TP=1$, $TP=2$,
and $TP=4$, all with $PP=1$. The left plot examines how average time spent
swapping scales with TP. These trials confirm that the swap latency does
decrease as TP increases, as we hypothesized. However, the latency on a single
GPU is noticeably higher than the lower bound, and the scaling appears to be
less than linear; this difference may be explained by the alpha-beta
communication model. Model parameters are transmitted not as one long stream,
but as separate messages for the individual tensors. Each TP shard still
contains the same number of tensors as the original model albeit smaller, so
the same number of messages must be sent by each worker when loading. Taking
the expression $\alpha+\beta*n$ for total latency, while message size $n$ is
reduced per worker, the per message latency term $\alpha$ remains the same,
leading to sublinear scaling.

The right plot shows swapping and execution times in proportion to the end-to-
end latency. From the plot, it is clear that swapping latency remains the
bottleneck in all cases, but as the number of GPUs increases, the proportion
of overall time spent swapping decreases; this highlights how model
parallelism benefits swapping even more than execution.

Figure 6: Swapping latency with changing PP scale.

Next, we run an experiment that varies PP degree between 1, 2, and 4 worker
stages. Similar to the TP experiments, Fig. 6 shows that increasing PP also
decreases the swapping latency. We postulate that in this case, sublinear
scaling stems from delays as a load entry is pipelined through workers. Since
workers process batch entries synchronously, load entries, despite being
asynchronous, must still wait for their turn.

Figure 7: Swapping latency for $TP=2$, $PP=2$.

TP and PP are often used together in practice, so we ran an additional trial
with $TP=2$, $PP=2$. From Fig. 7, we see that the mixed parallelism
configuration has lower latency than both pure TP and pure PP for the same
number of workers, and it in fact approaches the ideal scaling target. The
positive effect of mixing parallelism may be because the previously described
overheads in the TP case and in the PP case are lessened at smaller degrees.

###  5.2 Simulated Workloads

We characterize the practical performance of our model parallel swapping
design with simulated workloads of serving multiple OPT-13B models on the same
cluster of four A100 GPUs. For all experiments conducted, we follow the
configuration of $TP=2$, $PP=2$.

Each simulation trial begins with several warm up requests that are not
recorded. Then, requests are sent to all models over a 30 second period, with
the distribution of requests to each model following a random independent
Gamma arrival process. Each request has an input token length of 8. Across
simulations, we vary two parameters: the assignment of mean arrival rates to
each model, and the coefficient of variation (CV) that is shared by all
models. For our purposes, assigning different mean arrival rates simulates how
request rates skew toward a subset of models, while CV adjusts the burstiness
of requests. For instance, $CV=4$ is a high degree of burstiness, and
$Rates=(10,1,1)$ represents a skew toward the first model relative to the
other two.

Our first set of simulations serves three models at once, limiting to at most
two models in GPU memory at all times, and we check that GPU memory usage
approximately matches the footprint of two OPT-13B models. The maximum batch
size is 8. Tab. 1 summarizes the average end-to-end latencies for the grid of
parameters we measured, with three variations in skew and three variations in
CV. Fig. 8 visualizes request latency CDFs of all models combined for each
pair of (Rates, CV).

We observe a common pattern that as CV increases from 0.25 to 4, the latency
tends to decrease, which can be seen both in the table of average latencies
and in the CDF curves in each plot shifting toward the top left corner. This
suggests Computron performs better when request patterns are bursty.
Intuitively, bursty request distributions mean higher likelihood of
consecutive requests to the same model, so less model swaps occur because the
engine schedules the oldest request first with LRU replacement.

In the three-model simulation, changing the skew of request rates only
marginally increases the maximum latency, and in general has little impact on
the overall latency distribution. This provides evidence that Computron can
tolerate workloads with imbalanced requests rates. Though our simulations only
test static request rate assignments, we expect that this tolerance can also
extend to dynamic scenarios where the skewness of request rates changes over
time.

Skew | $CV=0.25$ | $CV=1$ | $CV=4$
---|---|---|---
$(1,1,1)$ | 1.262 | 0.606 | 0.518
$(10,1,1)$ | 1.172 | 0.886 | 0.550
$(10,10,1)$ | 1.014 | 0.716 | 0.374
Table 1: Average latencies for combinations of (Rates, CV) when serving 3
models with only 2 models in GPU memory. Figure 8: Latency CDFs for
combinations of (Rates, CV) when serving 3 models with only 2 models in GPU
memory.

The second set of simulations serves six models at once, limiting to at most
four models in GPU memory at all times and with the maximum batch size set to
32. The results in Fig. 9 show similar patterns as the previous three-model
simulations. When $CV=4$, the latency distribution of serving six models is
actually lower than serving three models on average based on Tab. 2, which
indicates that good resource utilization can be achieved when requests are
bursty. On the other hand, latencies of lower CV trials are scaled by
approximately a factor of two. A possible explanation for this is that in
lower CV trials, GPUs have already been maximally utilized conditioned on the
request distribution and scheduling order, so doubling the workload leads to
doubled latency.

Looking at the bigger picture, both the three-model and the six-model
simulations reveal that many requests take longer than the isolated latency
measurements from §5.1 would suggest. Two possible causes for this are the
scheduling algorithm and the choice of maximum batch size. Our simple oldest-
request-first scheduling algorithm overlooks global information that may be
used to reduce average latency. The maximum batch size trades off between the
rate at which a model’s request queue is drained and the compute time of that
batch, and it may have some interactions with the request arrival
distributions. We defer more thorough investigation of these effects to future
work.

Skew | $CV=0.25$ | $CV=1$ | $CV=4$
---|---|---|---
$(1,1,1,1,1,1)$ | 1.847 | 1.282 | 0.174
$(10,10,1,1,1,1)$ | 2.017 | 1.413 | 0.229
$(10,10,10,10,1,1)$ | 1.535 | 1.470 | 0.312
Table 2: Average latencies for combinations of (Rates, CV) when serving 6
models with only 4 models in GPU memory. Figure 9: Latency CDFs for
combinations of (Rates, CV) when serving 6 models with only 4 models in GPU
memory.

##  6 Conclusion

We design and implement Computron, a system that is capable of serving
multiple deep learning models with billions of parameters on the same cluster
of GPUs, and it can exceed aggregate GPU memory capacity through model
parallel swapping. In isolated tests, we demonstrate that our design takes
advantage of both TP and PP to speed up the swapping of distributed models.
Simulated random workloads show that our system can tolerate bursty and skewed
request patterns. These features enable organizations to efficiently serve
many cutting-edge large models for different tasks when compute resources are
limited.

An optimization that may significantly improve serving performance is to
speculatively load or offload models. In real world scenarios, requests to
different models are often not independent processes, but instead have
predictable patterns, such as the same model being requested many times
consecutively to generate a sequence, a subset of models often being requested
in some fixed order, or a model being more frequently requested at a certain
time of day. More sophisticated load scheduling algorithms with predictive
capabilities can drastically reduce the number of on-demand swaps, and by
extension, serving latency.

A problem that has not been resolved in this work is handling models with
different sizes, and even different model parallelism configurations. Our
system currently assumes that every model instance is evenly distributed
across the cluster in the same way and with the same memory footprint.
Removing that assumption brings many complexities such as the decision problem
of what to load/offload when swapping and whether workers should handle each
model differently.

## References

  * [1] Reza Yazdani Aminabadi, Samyam Rajbhandari, Ammar Ahmad Awan, Cheng Li, Du Li, Elton Zheng, Olatunji Ruwase, Shaden Smith, Minjia Zhang, Jeff Rasley, and Yuxiong He.  Deepspeed-inference: Enabling efficient inference of transformer models at unprecedented scale.  In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC ’22. IEEE Press, 2022\.
  * [2] Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language models are few-shot learners, 2020.
  * [3] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding, 2019.
  * [4] Jiangsu Du, Ziming Liu, Jiarui Fang, Shenggui Li, Yongbin Li, Yutong Lu, and Yang You.  Energonai: An inference system for 10-100 billion parameter transformer models, 2022.
  * [5] Arpan Gujarati, Reza Karimi, Safya Alzayat, Wei Hao, Antoine Kaufmann, Ymir Vigfusson, and Jonathan Mace.  Serving DNNs like clockwork: Performance predictability from the bottom up.  In 14th USENIX Symposium on Operating Systems Design and Implementation (OSDI 20), pages 443–462. USENIX Association, November 2020.
  * [6] Shenggui Li, Jiarui Fang, Zhengda Bian, Hongxin Liu, Yuliang Liu, Haichen Huang, Boxiang Wang, and Yang You.  Colossal-ai: A unified deep learning system for large-scale parallel training, 2022.
  * [7] Zhuohan Li, Lianmin Zheng, Yinmin Zhong, Vincent Liu, Ying Sheng, Xin Jin, Yanping Huang, Zhifeng Chen, Hao Zhang, Joseph E. Gonzalez, and Ion Stoica.  Alpaserve: Statistical multiplexing with model parallelism for deep learning serving, 2023.
  * [8] Deepak Narayanan, Mohammad Shoeybi, Jared Casper, Patrick LeGresley, Mostofa Patwary, Vijay Anand Korthikanti, Dmitri Vainbrand, Prethvi Kashinkunti, Julie Bernauer, Bryan Catanzaro, Amar Phanishayee, and Matei Zaharia.  Efficient large-scale language model training on gpu clusters using megatron-lm, 2021.
  * [9] Nersc.  Perlmutter architecture.  https://docs.nersc.gov/systems/perlmutter/architecture/.
  * [10] Long Ouyang, Jeff Wu, Xu Jiang, Diogo Almeida, Carroll L. Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, John Schulman, Jacob Hilton, Fraser Kelton, Luke Miller, Maddie Simens, Amanda Askell, Peter Welinder, Paul Christiano, Jan Leike, and Ryan Lowe.  Training language models to follow instructions with human feedback, 2022\.
  * [11] Shaden Smith, Mostofa Patwary, Brandon Norick, Patrick LeGresley, Samyam Rajbhandari, Jared Casper, Zhun Liu, Shrimai Prabhumoye, George Zerveas, Vijay Korthikanti, Elton Zhang, Rewon Child, Reza Yazdani Aminabadi, Julie Bernauer, Xia Song, Mohammad Shoeybi, Yuxiong He, Michael Houston, Saurabh Tiwary, and Bryan Catanzaro.  Using deepspeed and megatron to train megatron-turing nlg 530b, a large-scale generative language model, 2022.
  * [12] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, Aurelien Rodriguez, Armand Joulin, Edouard Grave, and Guillaume Lample.  Llama: Open and efficient foundation language models, 2023.
  * [13] Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, Todor Mihaylov, Myle Ott, Sam Shleifer, Kurt Shuster, Daniel Simig, Punit Singh Koura, Anjali Sridhar, Tianlu Wang, and Luke Zettlemoyer.  Opt: Open pre-trained transformer language models, 2022.
  * [14] Lianmin Zheng, Wei-Lin Chiang, Ying Sheng, Siyuan Zhuang, Zhanghao Wu, Yonghao Zhuang, Zi Lin, Zhuohan Li, Dacheng Li, Eric. P Xing, Hao Zhang, Joseph E. Gonzalez, and Ion Stoica.  Judging llm-as-a-judge with mt-bench and chatbot arena, 2023.


Quantification of Ebola virus replication kinetics in vitro

Laura E. Liaoa,1 Jonathan Carruthersb,1 Sophie J. Smitherc CL4 Virology
Teamc,2 Simon A. Wellerc Diane Williamsonc Thomas R. Lawsc Isabel García-
Dorivald Julian Hiscoxd Benjamin P. Holdere Catherine A. A. Beaucheminf,g Alan
S. Perelsona Martín López-Garcíab Grant Lytheb John Barrh Carmen Molina-
Parísb,3

a Theoretical Biology and Biophysics, Los Alamos National Laboratory, Los
Alamos, NM, USA 87545

b Department of Applied Mathematics, School of Mathematics, University of
Leeds, Leeds LS2 9JT, UK

c Defence Science and Technology Laboratory, Salisbury SP4 0JQ, UK

d Institute of Infection and Global Health, University of Liverpool,
Liverpool, L69 7BE, UK

e Department of Physics, Grand Valley State University, Allendale, MI, USA
49401

f Department of Physics, Ryerson University, Toronto, ON, Canada M5B 2K3

g Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS) Research
Program at RIKEN, Wako, Saitama, Japan, 351-0198

h School of Molecular and Cellular Biology, University of Leeds, Leeds LS2
9JT, UK

1 These authors contributed equally to this work.

2 Membership list can be found in the Acknowledgments section.

3<EMAIL_ADDRESS>
## Abstract

Mathematical modelling has successfully been used to provide quantitative
descriptions of many viral infections, but for the Ebola virus, which requires
biosafety level 4 facilities for experimentation, modelling can play a crucial
role. Ebola modelling efforts have primarily focused on in vivo virus
kinetics, e.g., in animal models, to aid the development of antivirals and
vaccines. But, thus far, these studies have not yielded a detailed
specification of the infection cycle, which could provide a foundational
description of the virus kinetics and thus a deeper understanding of their
clinical manifestation. Here, we obtain a diverse experimental data set of the
Ebola infection in vitro, and then make use of Bayesian inference methods to
fully identify parameters in a mathematical model of the infection. Our
results provide insights into the distribution of time an infected cell spends
in the eclipse phase (the period between infection and the start of virus
production), as well as the rate at which infectious virions lose infectivity.
We suggest how these results can be used in future models to describe co-
infection with defective interfering particles, which are an emerging
alternative therapeutic.

## Author summary

The two deadliest Ebola epidemics have both occurred in the past five years,
with one of these epidemics still ongoing. Mathematical modelling has already
provided insights into the spread of disease at the population level as well
as the effect of antiviral therapy in Ebola-infected animals. However, a
quantitative description of the replication cycle is still missing. Here, we
report results from a set of in vitro experiments involving infection with the
Ecran strain of Ebola virus. By parameterizing a mathematical model, we are
able to determine robust estimates for the duration of the replication cycle,
the infectious burst size, and the viral clearance rate.

## Introduction

The world’s second largest Ebola outbreak is currently underway in the
Democratic Republic of Congo. Ebola virus (EBOV) causes severe and fatal
disease with death rates of up to 90% [1]. There is an urgent need to prevent
and treat EBOV infections, but no antiviral drugs or monoclonal antibodies
have been approved in Africa, the EU, or the US. Recently the first EBOV
vaccine has been approved by European regulators [2]. Experimental therapies
[3], including antiviral drugs (remdesivir [4] and favipiravir [5, 6]) and a
cocktail of monoclonal antibodies (ZMapp) [7], have been assessed in the
2013–2016 West Africa Ebola virus disease outbreak. Other promising monoclonal
antibody therapies, called mAb114 and REGN-EB3, have been deployed in the
current 2018–2019 Kivu Ebola epidemic [8]. A better understanding of the
precise infection kinetics of EBOV is warranted.

Mathematical modelling of viral dynamics has provided a quantitative
understanding of within-host viral infections, such as HIV [9], influenza
[10], Zika [11], and more recently, EBOV. Mathematical modelling studies have
analyzed the plasma viral load dynamics of EBOV-infected animals (mice [12],
non-human primates [13, 14]) while under therapy with favipiravir, and have
identified estimates of favipiravir efficacy and target drug concentrations.
In addition, mechanistic models of innate and adaptive immune responses were
used to provide an explanation of EBOV infection dynamics in non-human
primates [14], and of differences between fatal and non-fatal cases of human
infection [15]. Moreover, mathematical models have been used to predict the
effect of treatment initiation time on indicators of disease severity [12, 15]
and survival rates [14], to predict the clearance of EBOV from seminal fluid
of survivors [16], and to theoretically explore treatment of EBOV-infected
humans with antivirals that possess different mechanisms of action (i.e.,
nucleoside analog, siRNA, antibody) [15].

Alongside the progress made in understanding within-host infections, a
complementary view of infection can be provided by mathematical modelling of
infections at the in vitro level. Combined with in vitro time course data,
mathematical models (MMs) have provided a detailed quantitative description of
the viral replication cycle of influenza A virus [17, 18], SHIV [19, 20], HIV
[21], and other viruses [22, 23, 24, 25, 26]. Such studies yield estimates of
key quantities such as the basic reproductive number (defined as the number of
secondary infections caused by one infected cell in a population of fully
susceptible cells), half-life of infected cells, and viral burst size, which
cannot be obtained directly from data [27]. In the context of in vitro
infections, parameterized MMs have been used to predict the outcome of
competition experiments between virus strains [28, 29, 30] (i.e., which strain
dominates in a mixed infection), map differences in genotype to changes in
phenotype [28, 30] (e.g., associate a single mutation to ten-fold faster viral
production), quantify fitness differences between virus strains [31, 32]
(e.g., which strain has a larger infectious burst size), quantify the
contribution of different modes of transmission (cell-to-cell versus cell-
free) [21], and identify the target of antiviral candidates [33] (e.g.,
whether a drug inhibits viral entry or viral production). One prior study [34]
utilized in vitro infection data from the literature to estimate EBOV
infection parameters, but had several parameter identifiability issues due to
insufficient data.

Our goal is to obtain robust estimates of viral infection parameters that
characterize the EBOV replication cycle. We follow a mathematical modelling
approach that has been successfully applied in the analysis of other viral
infections in vitro [17]. To this end, we performed a suite of in vitro
infection assays (single-cycle, multiple-cycle, and viral infectivity decay
assays) using EBOV and Vero cells, and collected detailed extracellular
infectious and total virus time courses. The viral kinetic data were simulated
with a multicompartment ordinary differential equation MM, and posterior
distributions of the MM parameters were estimated using a Markov chain Monte
Carlo (MCMC) approach. We estimate that one EBOV-infected cell spends
$\sim$30\text{\,}\mathrm{h}$$ in an eclipse phase before it releases
infectious virions at a rate of $13\text{\,}{\mathrm{h}}^{-1}$, over its
infectious lifetime of $\sim$83\text{\,}\mathrm{h}$$. The number of infectious
virions produced over an infected cell’s lifetime is $\sim 1000$, with an
estimated basic reproductive number of $\sim 600$. We also discuss challenges
in collecting other types of virus dynamic data (e.g., intracellular viral RNA
or cell counts).

## Results

#### Ebola Virus Kinetics In Vitro

Vero cell monolayers were infected with EBOV at a multiplicity of infection
(MOI) of 5, 1, . Infectious and total virus concentrations were determined
from extracellular virus harvested from the supernatant of each well at
various times post-infection (Fig. 1 A–C, E—F). At the start of infection, the
virus concentrations do not rise for some time, reflecting the time it takes
for viral entry, replication and release. After $24\text{\,}\mathrm{h}$, the
virus concentrations grow exponentially as infected cells begin producing
virus. When all cells in the well are infected, the virus concentrations peak
at approximately and
$10^{13}\text{\,}\mathrm{c}\mathrm{o}\mathrm{p}\mathrm{y}\mathrm{/}\mathrm{m}\mathrm{L}$
and the peak is sustained for $\sim$72\text{\,}\mathrm{h}$$. Thereafter, the
virus concentrations decline when virus production ceases, presumably due to
the death of infected cells. Additionally, the kinetics of viral infectivity
decay and virus degradation were assessed with a mock yield assay (Fig. 1 D,
G). In the mock yield assay, an inoculum of virus was incubated in wells under
the same conditions as the growth assays, but in the absence of cells, and
sampled over time.

Fig 1: Kinetics of EBOV infection in vitro and mock yield assays. Vero cell
monolayers were infected with EBOV at a multiplicity of infection (MOI) 5, 1,
or , as indicated. At various times post-infection, the infectious
($\text{TCID}_{50}$/mL; A–C) and total virus (copy/mL; E–G) in the supernatant
were determined. A mock yield assay was also performed to quantify the decay
of infectious (D) and total virus (G). In each assay, the experimental data
(circles) were collected either in duplicate (MOI 5) or triplicate (all other
assays). Note that the total virus concentration collected in the MOI 5
infection was omitted from the analysis due to inconsistencies in the peak
value (S1 Appendix, Fig. A). The lines represent the pointwise median of the
time courses simulated from our MM, which are bracketed by 68% (light grey)
and 95% (dark grey) credible regions (CR). These data were used to extract the
posterior probability likelihood distributions of the infection parameters
(Fig. 2). Note that parameters of the calibration curve used to convert cycle
threshold values (Ct) to total virus (copy/mL) were also estimated (Fig. 3).
The variability introduced from this conversion is shown by two error bars on
each total virus data point, indicating the 68% (same colour) and 95% (black)
CR.

#### Mathematical Model of Viral Infection and Parameter Estimates

The in vitro EBOV infection kinetics were captured with a MM that has been
used successfully in past works to capture influenza A virus infection
kinetics in vitro [30, 28, 31]. The MM is given by the system of ordinary
differential equations:

$\displaystyle\frac{\mathrm{d}T}{\mathrm{d}t}$ $\displaystyle=-\beta
TV_{\text{inf}}$ $\displaystyle\frac{\mathrm{d}E_{1}}{\mathrm{d}t}$
$\displaystyle=\beta TV_{\text{inf}}-\frac{n_{E}}{\tau_{E}}E_{1}$
$\displaystyle\frac{\mathrm{d}E_{i=2,3,...,n_{E}}}{\mathrm{d}t}$
$\displaystyle=\frac{n_{E}}{\tau_{E}}\left(E_{i-1}-E_{i}\right)$
$\displaystyle\frac{\mathrm{d}I_{1}}{\mathrm{d}t}$
$\displaystyle=\frac{n_{E}}{\tau_{E}}E_{n_{E}}-\frac{n_{I}}{\tau_{I}}I_{1}$
(1) $\displaystyle\frac{\mathrm{d}I_{j=2,3,...,n_{I}}}{\mathrm{d}t}$
$\displaystyle=\frac{n_{I}}{\tau_{I}}\left(I_{j-1}-I_{j}\right)$
$\displaystyle\frac{\mathrm{d}V_{\text{inf}}}{\mathrm{d}t}$
$\displaystyle=p_{\text{inf}}\sum_{j=1}^{n_{I}}I_{j}-c_{\text{inf}}V_{\text{inf}}$
$\displaystyle\frac{\mathrm{d}V_{\text{tot}}}{\mathrm{d}t}$
$\displaystyle=p_{\text{tot}}\sum_{j=1}^{n_{I}}I_{j}-c_{\text{tot}}V_{\text{tot}}$

In this MM, susceptible uninfected target cells $T$ can be infected by
infectious virus $V_{\text{inf}}$ with infection rate constant $\beta$, and
subsequently enter the non-productive eclipse phase $E_{i=1,\ldots,n_{E}}$,
followed by a transition into the productively infectious phase
$I_{j=1,\ldots,n_{I}}$. The eclipse and infectious phases are divided into a
number of compartments given by $n_{E}$ and $n_{I}$, respectively, such that
the time spent in each phase follows an Erlang distribution with an average
duration of $\tau_{E,I}\pm\frac{\tau_{E,I}}{\sqrt{n_{E,I}}}$. While cells are
in the infectious phase, they produce infectious (total) virus
$V_{\text{inf}}$ ($V_{\text{tot}}$) at a rate $p_{\text{inf}}$
($p_{\text{tot}}$), which lose infectivity (viability) at rate
$c_{\text{inf}}$ ($c_{\text{tot}}$). The MM Eq. (1) captures both infectious
virus, quantified by $\text{TCID}_{50}$ measurements of supernatant samples,
and total virus, quantified by quantitative, real-time, reverse transcriptase
PCR (hereafter, RT-qPCR). The latter experimental quantity was obtained by
converting cycle threshold (Ct) values from RT-qPCR to copy number (Fig. 3)
using Eq. (3) (Methods).

Predicted virus time courses from the MM are shown in Fig. 1 where the solid
lines represent the pointwise median and the grey bands show narrow 95%
credible regions (CR), indicating that the MM reproduces the viral kinetic
data well. Using a Markov chain Monte Carlo (MCMC) approach, we obtained
posterior probability likelihood distributions (PostPLDs) for each of the MM
parameters (Fig. 2). Narrow PostPLDs were extracted with mild correlations
between parameters (S1 Appendix, Fig. C), indicating good practical
identification of all parameters.

#### A Quantitative Description of the EBOV Lifecycle

The MCMC analysis gives us the following quantitative description of the EBOV
lifecycle within Vero cells. An EBOV-infected Vero cell spends approximately
$30\text{\,}\mathrm{h}$ with a 95% credible region of
[$26\text{\,}\mathrm{h}$, $37\text{\,}\mathrm{h}$] in the eclipse phase before
progeny EBOV successfully bud. Subsequently, infectious virus is produced at a
rate of 13 [10, 20] virions per cell per hour over a duration of
$83\text{\,}\mathrm{h}$ [$64\text{\,}\mathrm{h}$, $95\text{\,}\mathrm{h}$]
before virus production ceases due to cell death. An infectious burst size of
1096 progeny virions [1000, 1259] is released from each infected cell over its
virus-producing lifetime. Once infectious virus enters the cell culture
medium, they lose infectivity at a rate of $0.06\text{\,}{\mathrm{h}}^{-1}$
[$0.055\text{\,}{\mathrm{h}}^{-1}$, $0.068\text{\,}{\mathrm{h}}^{-1}$], which
is comparable to other viruses such as influenza A virus [31] or SHIV [19].
Overall, the in vitro spread of infection is rapid, as characterized by an
infecting time of $2\text{\,}\mathrm{h}$ [$1.6\text{\,}\mathrm{h}$,
$2.7\text{\,}\mathrm{h}$], which is defined as the time required for a single
virus-producing cell to infect one more [35]. These dynamics imply a large
basic reproductive number of $589\text{\,}\mathrm{[}$398, 1000], which is
defined as the number of secondary infections caused by a single infected cell
in a population of fully susceptible cells.

Notably, we find that the durations of both the eclipse and infectious phases
follow a normal-like distribution, as given by $n_{E}$ of 13 [8, 23] and
$n_{I}$ of 14 [3, 85]. This implies that the eclipse phase comprises a
sequence of many distinct steps of short duration, without any one step
lasting significantly longer than the rest. Likewise, the same interpretation
applies to the infectious phase. The normal-like distribution of the eclipse
phase resembles that of influenza A virus [30], but contrasts with the fat-
tailed eclipse phase distribution of SHIV [20] which is likely due to a
process in the phase that is longer than the rest (e.g., integration).
Moreover, neither the eclipse nor infectious phase are exponentially
distributed ($n=1$) as is commonly assumed in analyses with MMs. Such an
assumption has been shown to impact estimates of antiviral efficacy that are
based on patterns of viral load decay under simulated therapy in HIV patients
[20].

Fig 2: Estimated parameter distributions of EBOV infection in vitro. Posterior
probability likelihood distributions (PostPLDs) of parameters in the MM (A–G)
were estimated using MCMC and the data in Fig. 1. Secondary parameters were
derived from these estimates (H–J). Note that the PostPLDs corresponding to
the number of eclipse and infectious phase compartments are integer-valued.
The remaining PostPLDs of parameters describing the total virus and
calibration curve are in S1 Appendix, Fig. B.

## Discussion

In this work we performed time-course Ebola virus (EBOV) infection experiments
at multiple MOIs in vitro and applied MCMC methods to precisely parameterize a
mathematical model (MM) of the infection. We extracted fundamental quantities
concerning the timing and viral production of EBOV replication. This
theoretical-experimental approach maximized the output of the costly and
difficult experiments, which must be performed in biosafety level 4
facilities. Previous studies of the EBOV lifecycle rely on safer virus-like
particles [36]. The only previously known MM of EBOV infection in vitro [34]
is restricted in its use due to problems with parameter identifiability;
specifically, the existence of strong correlations between parameters, such as
the rates of virus degradation and virus production. By obtaining a more
complete set of experimental observations, we have provided the first detailed
quantitative characterization of EBOV infection kinetics.

Some of our estimates of timescales in the EBOV infection kinetics fill gaps
in the knowledge of this virus, while others expose some tension with prior
mathematical modelling work. The eclipse phase, excluded from the previous in
vitro MM [34], has been found to be a significant part of the replication
cycle. Lasting approximately $30\text{\,}\mathrm{h}$, it is longer than the
eclipse phase for influenza A virus and HIV infections in humans
(4–$24\text{\,}\mathrm{h}$) [10, 37]. Although the eclipse phase is included
in existing MMs of EBOV-infected animals, its duration has never been
estimated, and the assumed values used in these studies were considerably
shorter than the value we identify here [12, 14]. Moreover, the observation
that the length of the eclipse phase follows an Erlang distribution is
contrary to these previous MMs, where it has been represented more simply as
an exponentially distributed time. These MMs also fix the value of the decay
rate of infectious virus to ensure that other parameters remain identifiable
[12]. Here, the robust estimate of this decay rate demonstrates the benefit of
performing a mock yield assay. Existing MMs of in vivo EBOV infection in
humans and non-human primates provide considerably shorter estimates of the
infection cycle (12.5–$15.3\text{\,}\mathrm{h}$) compared to the estimate of
$114\text{\,}\mathrm{h}$ ($\tau_{E}+\tau_{I}$) obtained here [12, 15]. Such a
difference is likely attributed to the inclusion of an implicit immune
response in these in vivo models, thereby accounting for the enhanced
clearance of infected cells by immune cells, such as CD8+ T cells [38]. This
also explains why a faster viral decay rate can be expected in vivo, and
subsequently why estimates of the basic reproduction number are greater here
than those obtained from in vivo MMs (5.96-9.01) [12, 15]. It remains to be
determined whether Vero cells are representative of the cells targeted by EBOV
in vivo, but by understanding EBOV replication in Vero cells, we have a
foundation from which more complex cell culture models might be developed.

In addition to virus measurements, previous studies have included susceptible
and infected cell measurements to fully parameterize the MM and obtain robust
estimates of the viral kinetics parameters [19, 39]. We initially set out to
obtain a more diverse data set that also included the kinetics of dead cells
and intracellular RNA over the course of infection, but encountered unexpected
challenges. To quantify cell viability, we treated infected monolayers at
various times with Trypan blue, which stains cells that have lost the ability
to exclude dye. Unfortunately, we were unable to associate this marker of cell
death to a stage of the viral lifecycle in our MM without making additional
assumptions. Ultimately, when we extended the MM to include these data, the
newly introduced parameters were dependent on these assumptions, and the
extracted values of the original parameters were largely unaffected (S1
Appendix). To determine intracellular viral kinetics, the supernatants from
infected cell cultures were removed and the remaining monolayers were washed
and trypsinized for quantification via $\text{TCID}_{50}$ assay and RT-qPCR.
These samples showed a high level of EBOV RNA and $\text{TCID}_{50}$ as early
as 4 hours post-infection, which remained at a constant level up to 1 day
post-infection, but rose thereafter (S1 Appendix). Additionally, the ratio of
RNA-to-$\text{TCID}_{50}$ resembled the ratio observed in the supernatant.
Thus, these measurements likely reflect the large amount of cell-associated
virions that remained after washing, effectively obscuring the intracellular
RNA signal.

While a highly-controlled in vitro system was necessary to achieve our precise
characterization of the EBOV infection kinetics, the applicability of these
results to a clinical situation is not immediately obvious, and represents a
serious limitation of the study. Nevertheless, our findings have some
relevance to understanding the EBOV infection in vivo. EBOV initially
replicates within macrophages and dendritic cells in subcutaneous and
submucosal compartments, but dissemination in the blood results in the
infection of multiple organs throughout the body [40]. Many different cell
types are infected with varying susceptibility to infection, as well as
varying levels of viral replication. While the infection unfolds, EBOV blocks
IFN production early on [6, 41]. In this sense, studying the infection of Vero
cells—which are IFN-deficient—narrowly models the infection of one type of
epithelial cell during the early stages of an EBOV infection in vivo.

Vero cells serve as a standard host cell for replication and are widely used
for testing antivirals in vitro [42], as well as in the development of viral
vaccines [43, 44, 45]. Mathematical modelling of EBOV infections in vitro
using Vero cells has relevance to such applications, particularly in the study
of emerging therapeutics. While we provided a quantitative depiction of
extracellular infection by EBOV as a valuable first step, we envisioned that
the MM could be extended to include intracellular viral RNA kinetics had the
appropriate data been collected. Such multiscale modelling approaches have
been used to provide insight into virus growth and also to the understanding
of direct-acting antivirals [46, 47, 48, 49]. We hope that these experiences
might help guide future efforts to obtain informative cell and intracellular
data.

As an alternative antiviral strategy, there has been renewed interest in
pursuing defective interfering particles (DIPs) [50] of highly pathogenic
viruses. A DIP is a viral particle that contains defective interfering RNA (DI
RNA), which can be a shortened version of the parent genome that renders a DIP
replication-incompetent on its own (because it may lack the gene for an
essential viral component such as viral polymerase), but also elicits virus-
interfering properties. Within a cell co-infected by both DIPs and virus, the
DI RNA has a replicative advantage over the full-length RNA and outcompetes it
to produce more DIPs than virus progeny, effectively reducing the infectious
virus yield. EBOV DI RNA has been observed [51] but much remains to be
understood. Like with any other antiviral, MMs can be used to determine the
efficacy and mechanism of action of candidate DI RNAs, and to explore the
impact of dose and timing [15]. In particular, our estimates of EBOV infection
kinetics parameters are directly applicable to future mathematical modelling
of the interactions between EBOV and EBOV DIPs in vitro. Our estimated EBOV
infection parameters may also describe certain aspects of EBOV DIP infection.
For example, since DIPs have the same viral proteins and capsid as virions,
they would infect cells with the same infection rate constant, $\beta$. Since
DIPs also piggyback on the virus’ replication cycle, we might expect the same
eclipse and infectious phase lengths ($\tau_{E},\tau_{I}$) in a DIP and virus
co-infected cell.

In summary, the MM described here characterizes the replication cycle of EBOV
in a quantitative manner that will be beneficial for those creating in vitro
models to aid the development of antivirals and vaccines. We have made use of
a valuable set of in vitro results, carefully considering the structure of the
MM in order to maximize the information we can extract from them.

## Materials and Methods

### Cells and Virus

Vero C1008 cells (ECACC Cat. No.85020206) were obtained from Culture
Collection, Public Health England, UK. Vero C1008 cells were maintained in
Dulbecco’s minimum essential media supplemented with 10% (v/v) foetal calf
serum, 1% (v/v) L-glutamine and 1% (v/v) penicillin/streptomycin (Sigma). For
experimental purposes, the foetal calf serum concentration was reduced to 2%
(v/v).

Ebola virus H. sapiens-tc/COD/1976/Yambuku-Ecran, hereafter referred to as
EBOV was used in all studies. This virus, previously known as EBOV “E718” [52]
was supplied by Public Health England. Passage 5 material was used to infect
Vero C1008 cells. Virus was harvested on day 5 post-inoculation and titrated
to produce a working stock at .

### Quantification of Virus

EBOV was titrated in 96-well plates using the endpoint fifty percent tissue
culture infectious dose ($\text{TCID}_{50}$) assay [53]. Briefly, virus was
ten-fold serially diluted in 96 well plates of Vero C1008 cells. After one
week of incubation at $37\text{\,}\mathrm{\SIUnitSymbolDegree}\Celsius$/5%
$\mathrm{CO}_{2}$, all wells were observed under the microscope and scored for
presence or absence of cytopathic effects. The 50% endpoint was then
calculated using the method of Reed & Muench [54]. RNA extractions were
performed using the QiAMP Viral RNA Mini Kit (Qiagen, UK). Two
$50\text{\,}\mathrm{\SIUnitSymbolMicro L}$ elutions were performed for each
sample to increase the volume available for RT-PCR.

The genetic material of EBOV was quantified using the RealStar® Filovirus
Screen RT-PCR Kit (Altona diagnostics, Country) following the instructions of
the manufacturer. This assay has been performed many times against a standard
curve of plasmid containing the L gene from EBOV. The number of genomes can be
estimated from the Ct values as described in Eq. (3). In this context, the
number of genomes might consist of incomplete negative sense RNA molecules
encoding this sequence of the L gene. However, we do not believe that these
will be common ($<5\%$) based upon observations made with next generation
sequencing (paper in preparation). MOI 5 experiments were analysed using a
BIORAD CFX Connect – Real Time System, while samples for the remaining MOIs
were analysed using a QuantStudio 7 Flex Real-Time PCR System. Signal from
control RNA was compared between experiments and machines and we found no
evidence of differences. The parameters of the calibration curve required to
convert Ct values to total virus used samples from the MOI 5 experiments.

### Infections

Twenty-four-well plates were seeded with Vero C1008 cells at
$10^{5}\text{\,}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}\mathrm{/}\mathrm{m}\mathrm{L}$.
EBOV was added at MOIs of either 5, 1, or 0.1. Vero cells were grown to 90%
confluence for all infections. The cell culture medium was not changed during
the experiment and all cultures reached confluence within
$24\text{\,}\mathrm{h}$ (S1 Appendix). At pre-determined intervals post-
infection samples were taken by aspiration of supernatant from wells. Samples
were stored at $-80\text{\,}\mathrm{\SIUnitSymbolDegree}\Celsius$ prior to
enumeration by TCID50 assay and RNA extraction for PCR. Note that the RNA from
the MOI 5 infection was omitted from further analysis due to inconsistencies
in the peak viral RNA, compared to the MOI 1 and 0.1 infections (S1 Appendix,
Fig. A). The viability of Vero cells in the absence of infection is not known
under these conditions, however, we have observed these cells for 168 h at 24
h intervals and observed only occasional cells that can be stained with the
viability stain Trypan blue.

### Mock yield or infectivity decay assay

EBOV was added to twenty-four-well plates at a final estimated density of
$5\times 10^{5}$ $\text{TCID}_{50}$. At pre-determined intervals post-
infection samples were taken by aspiration of supernatant from wells. Samples
were stored at $-80\text{\,}\mathrm{\SIUnitSymbolDegree}\Celsius$ prior to
enumeration by $\text{TCID}_{50}$ assay and RNA extraction for PCR.

### Construction of the standard RT-qPCR curve

The concentration of viral genome copies (copy/mL) in a standard sample $i$
($V_{\text{STD},i}$) and the number of doubling RT-qPCR cycles
($C_{t,\text{STD},i}$) required for this concentration of copies to reach an
arbitrarily fixed, chosen threshold concentration ($Q_{t}$), are linked by the
equation

$\displaystyle Q_{t}$ $\displaystyle=V_{\text{STD},i}\
(2\varepsilon)^{C_{t,\text{STD},i}}$ $\displaystyle\ln(V_{\text{STD},i})$
$\displaystyle=\underbrace{\ln(Q_{t})}_{y\text{-intercept}}-\underbrace{\ln(2\varepsilon)}_{\text{slope}}\
C_{t,\text{STD},i}$ (2)

where $\varepsilon$ is the efficacy of the RT-qPCR doubling, which should
ideally be equal to one (i.e., exactly doubles at each cycle) but can vary
about this value. In constructing the standard curve, we took five standard
samples ($V_{\text{STD},i=1...5}$) with known copy concentrations (via their
mass) and determined their corresponding $C_{t,\text{STD},i}$. These data are
shown in Fig. 3.

Fig 3: Standard RT-qPCR curve. Cycle threshold values (Ct) were converted to
total virus (copy/mL) using the above calibration curve, where the parameters
of the curve were estimated as a part of the analysis. The lines represent the
pointwise median of the time courses simulated from our MM, which are
bracketed by 68% (light grey) and 95% (dark grey) CR.

### Conversion of sample RT-qPCR $C_{t}$ values into $V_{\text{tot}}$

In quantifying the concentration of total virus, $V_{\text{tot}}$ (copy/mL),
in the extracellular virus samples collected from infection experiments, Eq.
(2) was used as follows

$\displaystyle\ln({V_{\text{tot}}}_{,i})=\ln(Q_{t})-\ln(2\varepsilon)\
C_{t,\text{sample},i}\equiv\mathcal{F}(C_{t,\text{sample},i})$ (3)

where ${V_{\text{tot}}}_{,i}$ is the concentration of copies in sample $i$,
given its RT-qPCR-determined $C_{t,\text{sample},i}$ value. Here, $\ln(Q_{t})$
and $\ln(2\varepsilon)$ are two parameters to be estimated as part of the MCMC
parameter estimation process described later in this section. As different
values for these two parameters are sampled in the MCMC process, the total
virus concentration data points vary. The variation in the conversion is
denoted by error bars on each total virus data point in Fig. 1 E–G.

### Mock-yield assay model

Loss of virus infectivity or loss of viral genome integrity over time
typically follows an exponential decay [22]. As such, the mock-yield (MY) or
infectivity decay assay can be captured via $V(t)=V_{0}\,\mathrm{e}^{-c\,t}$,
such that the experimental MY data are expected to follow

$\displaystyle\ln(V_{\text{inf}}(t))$
$\displaystyle=\ln(V^{\text{MY}}_{\text{inf},0})-c_{\text{inf}}\ t$ (4)
$\displaystyle\ln(V_{\text{tot}}(t))$
$\displaystyle=\ln(V^{\text{MY}}_{\text{tot},0})-c_{\text{tot}}\ t$ (5)

where $V_{\text{inf}}(t)$ and $V_{\text{tot}}(t)$ are the concentrations of
infectious ($\text{TCID}_{50}$/mL) and total (copy/mL) virus after an
incubation of time $t$ under the same conditions used during the infection
experiments, given the EBOV rate of loss of infectivity ($c_{\text{inf}}$) or
integrity ($c_{\text{tot}}$), and initial concentrations,
$V^{\text{MY}}_{\text{inf},0}$ and $V^{\text{MY}}_{\text{tot},0}$. These data
are shown in Fig. 1 D, G.

### Simulated infections and parameter estimation

In estimating the MM parameters, the following experimental data were
considered simultaneously: the RT-qPCR standardized curve (5 data points), the
MY assays (24 data points: 4 time points in triplicate for $C_{t}$ and
$V_{\text{inf}}$), and three infection assays at MOI of 5 (24 data points: 6
time points in duplicate for $C_{t}$ and $V_{\text{inf}}$), MOI of 1 (54 data
points: 9 time points in triplicate for $C_{t}$ and $V_{\text{inf}}$), and MOI
of 0.1 (53 data points: 9 time points in triplicate for $C_{t}$ and
$V_{\text{inf}}$, minus one contaminated sample in $C_{t}$).

Eq. (2) was used to capture the RT-qPCR standard curve, and its agreement with
the 5 experimental data points was computed as the sum-of-squared residuals
(SSR)

$\displaystyle\text{SSR}^{\text{STD}}=\frac{\sum_{i=1}^{5}\left[\mathcal{F}(C_{t,\text{STD},i})-\ln(V_{\text{STD},i})\right]^{2}}{\sigma_{V_{\text{tot}}}^{2}}\
,$

where $\sigma_{V_{\text{tot}}}^{2}$ is the variance, or squared of the
standard error, in experimentally measured $V_{\text{tot}}$, which will be
discussed in more details below. Eq. (4) and Eq. (5) were used to capture the
MY experiment, performed in triplicate, and sampled at 4 time points, for each
of $C_{t}$ and $V_{\text{inf}}$, and agreement was computed as

$\displaystyle\text{SSR}^{\text{MY}}=\frac{\sum_{i=1}^{12}\left[\ln(V^{\text{MY}}_{\text{inf},0})-c_{\text{inf}}\
t_{i}-\ln(V_{\text{inf}}(t_{i}))\right]^{2}}{\sigma_{V_{\text{inf}}}^{2}}$
$\displaystyle+\frac{\sum_{i=1}^{12}\left[\ln(V^{\text{MY}}_{\text{tot},0})-c_{\text{tot}}\
t_{i}-\mathcal{F}(C_{t,i})\right]^{2}}{\sigma_{V_{\text{tot}}}^{2}}$

Finally, MM Eq. (1) was used to reproduce the infection experiments, performed
in triplicate, and at 3 different MOIs (5, 1, and 0.1), and quantified via
both $\text{TCID}_{50}$ and RT-qPCR. In reproducing the infections, initial
conditions (at $t=0$) were such that $T(0)=1$,
$E_{i=1,...,n_{E}}(0)=I_{j=1,...,n_{I}}(0)=0$, and the initial infectious and
total virus concentrations for the 3 MOIs were computed as:

$\displaystyle V_{\text{inf}}(0)$
$\displaystyle=V^{\text{INF}}_{\text{inf},0}\times\text{MOI}$ $\displaystyle
V_{\text{tot}}(0)$
$\displaystyle=V^{\text{INF}}_{\text{tot},0}\times\text{MOI}$

where MOI was either 5, 1 or 0.1, and ($V^{\text{INF}}_{\text{inf},0}$,
$V^{\text{INF}}_{\text{tot},0}$) are 2 parameters to be estimated. Agreement
between MM Eq. (1) and experimental infection data was computed as

$\displaystyle\text{SSR}^{\text{INF}}=\sum_{\text{MOI}=[5,1,0.1]}\frac{\sum_{i}\left[\ln(V_{\text{inf}}^{\text{MM}}(t_{i})-\ln(V_{\text{inf}}(t_{i}))\right]^{2}}{\sigma_{V_{\text{inf}}}^{2}}$
$\displaystyle+\frac{\sum_{i}\left[\ln(V_{\text{tot}}^{\text{MM}}(t_{i})-\mathcal{F}(C_{t,i})\right]^{2}}{\sigma_{V_{\text{tot}}}^{2}}$

where $\sigma_{V_{\text{inf}}}^{2}=0.1$ and $\sigma_{V_{\text{tot}}}^{2}=0.1$
correspond to the variance in $\ln(V_{\text{inf}})$ and $\ln(V_{\text{tot}})$,
respectively, estimated as the variance of the residuals between the 2 to 3
replicates of $\ln(V_{\text{inf}})$ or $\ln(V_{\text{tot}})$ measured at each
time point and their corresponding mean, across all (STD, MY, and INF)
experimental data collected.

A total of 15 parameters — 6 parameters associated with experimental
conditions ($\ln(Q_{t})$, $\ln(2\varepsilon)$,
$\ln(V^{\text{MY}}_{\text{inf},0})$, $\ln(V^{\text{MY}}_{\text{tot},0})$,
$V^{\text{INF}}_{\text{inf},0}$, $V^{\text{INF}}_{\text{tot},0}$) and 9
parameters more closely associated with EBOV infection kinetics
($c_{\text{inf}}$, $c_{\text{tot}}$, $p_{\text{inf}}$, $p_{\text{tot}}$,
$\beta$, $\tau_{E}$, $\tau_{I}$, $n_{E}$, $n_{I}$) — were estimated (Table 1)
from 160 experimental data points using the python MCMC implementation phymcmc
[55], a wrapping library for emcee [56]. Posterior probability likelihood
distributions (PostPLDs) were obtained based on the parameter likelihood
function

$\displaystyle\ln(\mathcal{L}(\vec{p}))=-\frac{1}{2}\left[\text{SSR}^{\text{STD}}(\vec{p})+\text{SSR}^{\text{MY}}(\vec{p})+\text{SSR}^{\text{INF}}(\vec{p})\right]$

and the assumption of linearly uniform or $\ln$-uniform priors, where
$\vec{p}$ is the 15-parameter vector.

Table 1: Estimated parameters of EBOV infection in vitro. Parameter | Mode [95% CR]
---|---
Infectiousness, $\beta$ ($\frac{\mathrm{mL}}{\text{TCID}_{50}\cdot\mathrm{h}}$) | $10^{-6.48\ [-6.7,-6.3]}$
Eclipse phase length, $\tau_{E}$ (h) | $30.5\ [26,37]$
Number of eclipse compartments, $n_{E}$ | $13\ [8,23]$
Infectious phase length, $\tau_{I}$ (h) | $83.2\ [64,95]$
Number of infectious compartments, $n_{I}$ | $14\ [3,85]$
Infectious virus production rate, $p_{\text{inf}}$ ($\frac{\text{TCID}_{50}}{\mathrm{cell}\cdot\mathrm{h}}$) | $10^{1.12\ [1,1.3]}$
Total virus production rate, $p_{\text{tot}}$ ($\frac{\mathrm{RNA}}{\mathrm{cell}\cdot\mathrm{h}}$) | $10^{6.46\ [6.3,6.7]}$
Rate of loss of infectious virus, $c_{\text{inf}}$ (/h) | $0.0614\ [0.055,0.068]$
Rate of virus degradation, $c_{\text{tot}}$ (/h) | $0.00817\ [0.0035,0.013]$
Initial infectious virus inoculum, $V^{\text{INF}}_{\text{inf},0}$ ($\frac{\text{TCID}_{50}}{\mathrm{mL}}$) | $10^{5.39\ [5.3,5.5]}$
Initial total virus inoculum, $V^{\text{INF}}_{\text{tot},0}$ ($\frac{\mathrm{RNA}}{\mathrm{mL}}$) | $10^{11.5\ [11,12]}$
MY initial infectious virus inoculum, $\ln(V^{\text{MY}}_{\text{inf},0})$ | $13.7\ [13,14]$
MY initial total virus inoculum, $\ln(V^{\text{MY}}_{\text{tot},0})$ | $28.2\ [28,29]$
Standard RT-qPCR curve $y$-intercept, $\ln(Q_{t})$ | $37.8\ [37,39]$
Standard RT-qPCR curve slope, $\ln(2\varepsilon)$ | $0.613\ [0.57,0.66]$
Basic reproductive number, $R_{0}$ | $10^{2.77\ [2.6,3]}$
Infectious burst size, $p_{\mathrm{inf}}\tau_{I}$ ($\frac{\text{TCID}_{50}}{\mathrm{cell}}$) | $10^{3.04\ [3,3.1]}$
Infecting time, $t_{\mathrm{inf}}$ (h) | $10^{0.335\ [0.21,0.43]}$

## Supporting information

##### S1 Appendix.

Kinetics of cell-associated virus and cell viability.

## Acknowledgments

Members of the CL4 Virology Team include Lin Eastaugh, Lyn M. O’Brien, James
S. Findlay, Mark S. Lever, Amanda Phelps, Sarah Durley-White, Jackie Steward
and Ruth Thom. The authors would like to acknowledge Joseph Gillard for
helpful discussions and the International Centre for Mathematical Sciences
(ICMS), where the mathematical model was developed during a Research-in-Groups
programme.

## References

  * 1. Feldmann H, Geisbert TW.  Ebola haemorrhagic fever.  Lancet. 2011;377(9768):849–862.  doi:10.1016/S0140-6736(10)60667-8.
  * 2. Callaway E.  Make Ebola a thing of the past: First vaccine against deadly virus approved.  Nature. 2019;doi:10.1038/d41586-019-03490-8.
  * 3. Rojek A, Horby P, Dunning J.  Insights from clinical research completed during the west Africa Ebola virus disease epidemic.  Lancet Infect Dis. 2017;17(9):e280–e292.  doi:10.1016/S1473-3099(17)30234-7.
  * 4. Dörnemann J, Burzio C, Ronsse A, Sprecher A, De Clerck H, Van Herp M, et al.  First newborn baby to receive experimental therapies survives Ebola virus disease.  J Infect Dis. 2017;215(2):171–174.  doi:10.1093/infdis/jiw493.
  * 5. Sissoko D, Laouenan C, Folkesson E, M’Lebing AB, Beavogui AH, Baize S, et al.  Experimental treatment with favipiravir for Ebola virus disease (the JIKI trial): A historically controlled, single-arm proof-of-concept trial in Guinea.  PLoS Med. 2016;13(3):e1001967.  doi:10.1371/journal.pmed.1001967.
  * 6. Perez SC, Folkesson E, Anglaret X, Beavogui AH, Berbain E, Camara AM, et al.  Challenges in preparing and implementing a clinical trial at field level in an Ebola emergency: A cast study in Guinea, West Africa.  Plos Negl Trop Dis. 2017;11(6):e0005545.  doi:10.1371/journal.pntd.0005545.
  * 7. Group TPIW.  A randomized, controlled trial of ZMapp for Ebola virus infection.  N Engl J Med. 2016;375(15):1448–1456.  doi:10.1056/NEJMoa1604330.
  * 8. Maxmen A.  Experimental Ebola drugs face tough test in war zone.  Nature. 2018;561(7721):14.  doi:10.1038/d41586-018-06132-7.
  * 9. Perelson AS, Neumann AU, Markowitz M, Leonard JM, Ho DD.  HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time.  Science. 1996;271(5255):1582–1586.  doi:10.1126/science.271.5255.1582.
  * 10. Baccam P, Beauchemin CAA, Macken CA, Hayden FG, Perelson AS.  Kinetics of influenza A virus infection in humans.  J Virol. 2006;80(15):7590–7599.  doi:10.1128/JVI.01623-05.
  * 11. Best K, Guedj J, Madelain V, de Lamballerie X, Yon Lim S, Osuna CE, et al.  Zika plasma viral dynamics in nonhuman primates provides insights into early infection and antiviral strategies.  Proc Natl Acad Sci USA. 2017;114(33):8847–8852.  doi:10.1073/pnas.1704011114.
  * 12. Madelain V, Oestereich L, Graw F, Nguyen THT, de Lamballerie X, Mentré F, et al.  Ebola virus dynamics in mice treated with favipiravir.  Antiviral Res. 2015;123:70–77.  doi:10.1016/j.antiviral.2015.08.015.
  * 13. Guedj J, Piorkowski G, Jacquot F, Madelain V, Nguyen T, Rodallec A, et al.  Antiviral efficacy of favipiravir against Ebola virus: A translational study in cynomolgous macaques.  PLoS Med. 2018;15(3):e1002535.  doi:10.1371/journal.pmed.1002535.
  * 14. Madelain V, Baize S, Jacquot F, Reynard S, Fizet A, Barron S, et al.  Ebola viral dynamics in nonhuman primates provides insights into virus immuno-pathogenesis and antiviral strategies.  Nat Commun. 2018;9(1):4013.  doi:10.1038/s41467-018-06215-z.
  * 15. Martyushev A, Nakaoka S, Sato K, Noda T, Iwami S.  Modelling Ebola virus dynamics: Implications for therapy.  Antiviral Res. 2016;135:62–73.  doi:10.1016/j.antiviral.2016.10.004.
  * 16. Sissoko D, Duraffour S, Kerber R, Kolie JS, Beavogui AH, Camara AM, et al.  Persistence and clearance of Ebola virus RNA from seminal fluid of Ebola virus disease survivors: a longitudinal analysis and modelling study.  The Lancet Glob Health. 2017;5(1):e80–e88.  doi:10.1016/S2214-109X(16)30243-1.
  * 17. Handel A, Liao LE, Beauchemin CAA.  Progress and trends in mathematical modelling of influenza A virus infections.  Curr Opin Syst Biol. 2018;12:30–36.  doi:10.1016/j.coisb.2018.08.009.
  * 18. Möhler L, Flockerzi D, Sann H, Reichl U.  Mathematical model of influenza A virus production in large-scale microcarrier culture.  Biotechnol Bioeng. 2005;90(1):46–58.  doi:10.1002/bit.20363.
  * 19. Iwami S, Holder BP, Beauchemin CAA, Morita S, Tada T, Sato K, et al.  Quantification system for the viral dynamics of a highly pathogenic simian/human immunodeficiency virus based on an in vitro experiment and a mathematical model.  Retrovirology. 2012;9:18.  doi:10.1186/1742-4690-9-18.
  * 20. Beauchemin CAA, Miura T, Iwami S.  Duration of SHIV production by infected cells is not exponentially distributed: Implications for estimates of infection parameters and antiviral efficacy.  Sci Rep. 2017;7:42765.  doi:10.1038/srep42765.
  * 21. Iwami S, Takeuchi JS, Nakaoka S, Mammano F, Clavel F, Inaba H, et al.  Cell-to-cell infection by HIV contributes over half of virus infection.  Elife. 2015;4.  doi:10.7554/eLife.08150.
  * 22. Beauchemin CAA, Kim YI, Yu Q, Ciaramella G, DeVincenzo JP.  Uncovering critical properties of the human respiratory syncytial virus by combining in vitro assays and in silico analyses.  PLoS ONE. 2019;14(4):e0214708.  doi:10.1371/journal.pone.0214708.
  * 23. Gonzàlez-Parra G, De Ridder F, Huntjens D, Roymans D, Ispas G, Dobrovolny HM.  A comparison of RSV and influenza in vitro kinetic parameters reveals differences in infecting time.  PLoS ONE. 2018;13(2):e0192645.  doi:10.1371/journal.pone.0192645.
  * 24. Fukuhara M, Iwami S, Sato K, Nishimura Y, Shimizu H, Aihara K, et al.  Quantification of the dynamics of enterovirus 71 infection by experimental-mathematical investigation.  J Virol. 2013;87(1):701–705.  doi:10.1128/JVI.01453-12.
  * 25. Gonzàlez-Parra G, Dobrovolny HM, Aranda DF, Chen-Charpentier B, Rojas RAG.  Quantifying rotavirus kinetics in the REH tumor cell line using in vitro data.  Virus Res. 2018;244:53–63.  doi:10.1016/j.virusres.2017.09.023.
  * 26. Wethington D, Harder O, Uppulury K, Stewart WCL, Chen P, Kang T, et al.. Mathematical modeling identifies the role of adaptive immunity as a key controller of respiratory syncytial virus (RSV) titer in cotton rats; 2018.
  * 27. Iwami S, Sato K, Boer RJD, Aihara K, Miura T, Koyanagi Y.  Identifying viral parameters from in vitro cell cultures.  Front Microbiol. 2012;3:319.  doi:10.3389/fmicb.2012.00319.
  * 28. Pinilla LT, Holder BP, Abed Y, Boivin G, Beauchemin CAA.  The H275Y neuraminidase mutation of the pandemic A/H1N1 virus lengthens the eclipse phase and reduces viral output of infected cells, potentially compromising fitness in ferrets.  J Virol. 2012;86(19):10651–10660.  doi:10.1128/JVI.07244-11.
  * 29. Song H, Pavlicek JW, Cai F, Bhattacharya T, Li H, Iyer SS, et al.  Impact of immune escape mutations on HIV-1 fitness in the context of the cognate transmitted/founder genome.  Retrovirology. 2012;9:89.  doi:10.1186/1742-4690-9-89.
  * 30. Paradis EG, Pinilla L, Holder BP, Abed Y, Boivin G, Beauchemin CAA.  Impact of the H275Y and I223V mutations in the neuraminidase of the 2009 pandemic influenza virus in vitro and evaluating experimental reproducibility.  PLoS ONE. 2015;10(5):e0126115.  doi:10.1371/journal.pone.0126115.
  * 31. Simon PF, de La Vega MA, Paradis E, Mendoza E, Coombs KM, Kobasa D, et al.  Avian influenza viruses that cause highly virulent infections in humans exhibit distinct replicative properties in contrast to human H1N1 viruses.  Sci Rep. 2016;6:24154.  doi:10.1038/srep24154.
  * 32. Iwanami S, Kakizoe Y, Morita S, Miura T, Nakaoka S, Iwami S.  A highly pathogenic simian/human immunodeficiency virus effectively produces infectious virions compared with a less pathogenic virus in cell culture.  Theor Biol Med Model. 2017;14(1):9.  doi:10.1186/s12976-017-0055-8.
  * 33. Ikeda H, Godinho-Santos A, Rato S, Vanwalscappel B, Clavel F, Aihara K, et al.  Quantifying the antiviral effect of IFN on HIV-1 replication in cell culture.  Sci Rep. 2015;5:11761.  doi:10.1038/srep11761.
  * 34. Nguyen VK, Binder SC, Boianelli A, Meyer-Hermann M, Hernandez-Vargas EA.  Ebola virus infection modeling and identifiability problems.  Front Microbiol. 2015;6:257.  doi:10.3389/fmicb.2015.00257.
  * 35. Holder BP, Beauchemin CAA.  Exploring the effect of biological delays in kinetic models of influenza within a host or cell culture.  BMC Public Health. 2011;11 Suppl 1:S10.  doi:10.1186/1471-2458-11-S1-S10.
  * 36. Biedenkopf N, Hoenen T.  Modeling the Ebolavirus life cycle with transcription and replication-competent viruslike particle assays.  In: Ebolaviruses. Springer; 2017. p. 119–131.
  * 37. Dixit NM, Markowitz M, Ho DD, Perelson AS.  Estimates of intracellular delay and average drug efficacy from viral load data of HIV-infected individuals under antiretroviral therapy.  Antivir Ther. 2004;9(2):237–246.
  * 38. Gupta M, Greer P, Mahanty S, Shieh WJ, Zaki SR, Ahmed R, et al.  CD8-mediated protection against Ebola virus infection is perforin dependent.  The Journal of Immunology. 2005;174(7):4198–42020.
  * 39. Schulze-Horsel J, Schulze M, Agalaridis G, Genzel Y, Reichl U.  Infection dynamics and virus-induced apoptosis in cell culture-based influenza vaccine production-Flow cytometry and mathematical modeling.  Vaccine. 2009;27:2712–2722.  doi:10.1016/j.vaccine.2009.02.027.
  * 40. Chertow DS, Shekhtman L, Lurie Y, Davey RT, Heller T, Dahari H.  Modeling challenges of Ebola virus–host dynamics during infection and treatment.  Viruses. 2020;12(1):106.
  * 41. Edwards MR, Liu G, Mire CE, Sureshchandra S, Luthra P, Yen B, et al.  Differential regulation of interferon responses by Ebola and Marburg virus VP35 proteins.  Cell reports. 2016;14(7):1632–1640.
  * 42. Postnikova E, Cong Y, DeWald LE, Dyall J, Yu S, Hart BJ, et al.  Testing therapeutics in cell-based assays: Factors that influence the apparent potency of drugs.  PLoS ONE. 2018;13(3):e0194880.  doi:10.1371/journal.pone.0194880.
  * 43. Barrett NP, Mundt W, Kistner O, Howard MK.  Vero cell platform in vaccine production: moving towards cell culture-based viral vaccines.  Expert Rev Vaccines. 2009;8(5):607–618.  doi:10.1586/erv.09.19.
  * 44. Barrett NP, Terpening SJ, Snow D, Cobb RR, Kistner O.  Vero cell technology for rapid development of inactivated whole virus vaccines for emerging viral diseases.  Expert Rev Vaccines. 2017;16(9):883–894.
  * 45. Paillet C, Forno G, Kratje R, Etcheverrigaray M.  Suspension-Vero cell cultures as a platform for viral vaccine production.  Vaccine. 2009;27(46):6464–6467.
  * 46. Guedj J, Dahari H, Rong L, Sansone ND, Nettles RE, Cotler SJ, et al.  Modeling shows that the NS5A inhibitor daclatasvir has two modes of action and yields a shorter estimate of the hepatitis C virus half-life.  Proc Natl Acad Sci USA. 2013;110(10):3991–3996.  doi:10.1073/pnas.1203110110.
  * 47. Heldt FS, Frensing T, Pflugmacher A, Gröpler R, Peschel B, Reichl U.  Multiscale modeling of influenza A virus infection supports the development of direct-acting antivirals.  PLoS Comput Biol. 2013;9:e1003372.  doi:10.1371/journal.pcbi.1003372.
  * 48. de M Quintela B, Conway JM, Hyman JM, Guedj J, dos Santos RW, Lobosco M, et al.  A new age-structured multiscale model of the hepatitis C virus life-cycle during infection and therapy with direct-acting antiviral agents.  Front Microbiol. 2018;9:601.  doi:10.3389/fmicb.2018.00601.
  * 49. Zitzmann C, Kaderali L.  Mathematical analysis of viral replication dynamics and antiviral treatment strategies: From basic models to age-based multi-scale modeling.  Front Microbiol. 2018;9:1546.  doi:10.3389/fmicb.2018.01546.
  * 50. Rezelj VV, Levi LI, Vignuzzi M.  The defective component of viral populations.  Curr Opin Virol. 2018;33:74–80.  doi:10.1016/j.coviro.2018.07.014.
  * 51. Calain P, Monroe MC, Nichol ST.  Ebola virus defective interfering particles and persistent infection.  Virology. 1999;262(1):114–128.  doi:10.1006/viro.1999.9915.
  * 52. Kuhn JH, Lofts LL, Kugelman JR, Smither SJ, Lever MS, Groen Gv, et al.  Reidentification of Ebola virus E718 and ME as Ebola virus/H.sapiens-tc/COD/1976/Yambuku-Ecran.  Genome Announc. 2014;2(6):pii: e01178–14.  doi:10.1128/genomeA.01178-14.
  * 53. Smither SJ, Lear-Rooney C, Biggins J, Pettitt J, Lever MS, Jr GGO.  Comparison of the plaque assay and 50% tissue culture infectious dose assay as methods for measuring filovirus infectivity.  J Virol Methods. 2013;193(2):565–571.  doi:10.1016/j.jviromet.2013.05.015.
  * 54. Reed LJ, Muench H.  A simple method of estimating fifty per cent endpoints.  Am J Epidemiol. 1938;27(3):493–497.  doi:10.1093/oxfordjournals.aje.a118408.
  * 55. Beauchemin CAA. phymcmc: A convenient wrapper for emcee; 2019.  https://github.com/cbeauc/phymcmc.
  * 56. Foreman-Mackey D, Hogg DW, Lang D, Goodman J.  emcee: The MCMC hammer.  Publ Astron Soc Pac. 2013;125(925):306–312.  doi:10.1086/670067.