key: cord-0426567-9jvygedl
authors: Landry, Nicholas
title: Effect of time-dependent infectiousness on epidemic dynamics
date: 2021-06-18
journal: nan
DOI: 10.1103/physreve.104.064302
sha: 54b503a0ece654aff812821b8838b9bf1c5b1ad0
doc_id: 426567
cord_uid: 9jvygedl

In contrast to the common assumption in epidemic models that the rate of infection between individuals is constant, in reality, an individual's viral load determines their infectiousness. We compare the average and individual reproductive numbers and epidemic dynamics for a model incorporating time-dependent infectiousness and a standard SIR model for both fully-mixed and category-mixed populations. We find that the reproductive number only depends on the total infectious exposure and the largest eigenvalue of the mixing matrix and that these two effects are independent of each other. When we compare our time-dependent mean-field model to the SIR model with equivalent rates, the epidemic peak is advanced and modifying the infection rate function has a strong effect on the time dynamics of the epidemic. We also observe behavior akin to a traveling wave as individuals transition through infectious states.

Epidemic modeling has a rich tradition in network science [1] [2] [3] [4] with standard models such as the SIS (Susceptible -Infected -Susceptible) and SIR (Susceptible -Infected -Removed) models, for which rigorous mathematical theory has been developed. There are also more complex spatio-temporal models that more accurately capture the dynamics of disease spread in the real world [5] . Much interest has been devoted to the accurate prediction of the spread of the SARS-CoV-2 pandemic [6, 7] and to answering questions such as the efficacy of different prevention measures and the risk factors of different social situations [8] [9] [10] . In traditional literature, the SIR model is a canonical example of modeling the spread of disease with total immunity. This model has common extensions such as the SEIR (Susceptible -Exposed -Infected -Recovered) when one wants to incorporate a latent period which captures delays between transmission and infectiousness. With most of these models, however, a key assumption is that an individual's infectivity is constant. However, we know that an individual's infectiousness varies over the duration of the infection, according to their viral load [11, 12] . We define a framework to extend the SIR model by dividing the single infectious compartment into n stages as has been considered by Ref. [13] , known as the SI K R model in Ref. [14] , and assigning each stage a different infection rate as in Refs. [15, 16] . Other approaches have been considered, such as the message-passing approach [17, 18] , mapping an individual's viral load to an infection probability [10] , and looking at an infection density function [14, 19] . We use this approach to examine fullymixed populations and theoretical networks constructed from category-based mixing, both static and temporal.

The structure of the paper is as follows. In Section II we describe a framework for modeling time-dependent infectiousness. In Section III we use this model to create theoretical predictions for the reproductive number, * nicholas.landry@colorado.edu apply these predictions to several common cases, and validate our theory with numerical simulations. Lastly, in Section IV we discuss the implications of our theory.

We propose a general mean-field model to describe the spread of an epidemic including time-dependent infectiousness. In the following, we will refer to this model as the viral load (VL) model.

We consider a population of N nodes. We assume that a node i's intrinsic infectiousness is solely determined by the amount of time it has been infected, τ , and its corresponding viral load at that time, denoted v i (τ ), although other factors may be involved as well [8] . Several studies have examined the correspondence between an individual's viral load and their infectiousness [11, 12] but for this study, we simply define β i (τ ), the infectious rate function, as the rate at which node i transmits infection having been infected for a duration of time τ . Note that in the case where an infectious threshold exists [10, 20] , we can express the function as β i (τ )I τ ∈δ , where δ = {τ | β i (τ ) ≥ η} and η is the infectious threshold. This infectious rate function can vary in response to many factors such as asymptomatic versus symptomatic infection or severity of symptoms and be considered as being drawn according to some distribution. For much of this study, however, we assume that while β i (τ ) is heterogeneous in time, that every member of the population has the same infectious rate function, i.e., β i (τ ) = β(τ ), i = 1 . . . N , though we relax this assumption later. We assume that nodes start in the susceptible compartment (S) and that an infected individual infected for time τ infects a susceptible node with rate β(τ ). We approximate β(τ ) by evaluating it at n discrete times τ j = j∆τ , where ∆τ is fixed and n∆τ = τ R , the recovery time. Then we divide the infectious compartment, I into n stages, I j , j = 1 . . . n, each with an associated infection rate β j , in a similar manner to Refs. [13, 16] . Lastly, nodes that transition through all infection states accumulate in the recovered (R) compartment. We assume that the flow of infected individual between subsequent infectious compartments is deterministic and that upon entering the first infectious stage, an individual passes through all the subsequent stages as shown in Fig. 1 , meaning that γ i = 1/∆τ where ∆τ = τ R /n.

In the following, we define the m-th moment of a quantity q as q m = N i=1 q m i /N when q is a discrete quantity and as q m (τ ) =

There are many studies exploring the effect of more realistic infectious behavior. In Ref. [13] , the authors use the n-stage SI K R model with constant infectiousness on a fully-mixed network so that the infectious waiting time is gamma-distributed. In Ref. [16] , the authors explore the SI K R model with variable infectiousness for fully-mixed networks. For both of these models, the authors allow healing and recovery to occur at every infectious stage. In Ref. [21] , the authors explore the SI K R link-closure model with a constant infection rate and solely consider static networks. They simulate their model numerically on homogeneous and Erdös-Rényi networks. In Ref. [17] , the authors consider a message-passing approach to model time-dependent infectiousness and simulate their results on a static network. In Ref. [18] , the authors present a non-Markovian edge-based compartment model, prove its equivalence to the message-passing model, and describe how other models compare to the message-passing approach. In Refs. [13, 16] the authors solely consider the fully-mixed case and in Refs. [17, 18, 21] the authors solely consider static networks. In contrast, our approach encompasses fullymixed, static, and temporal networks. In Refs. [13, 21] , though the authors consider an SI K R model, they specify that the infectious rate is constant in contrast to our model where we allow the rate to vary over time. In addition, in Refs. [13, 16, 21] , they assume Markovian transitions between infectious states in contrast to our approach which enforces deterministic transitions between infectious states (as in Ref. [10] ).

We derive the reproductive number for the viral load model described above that has been cast as a system of mean-field ODEs. First, we derive the reproductive number for a fully-mixed model and second, we derive the reproductive number for an arbitrary category-mixed population. We comment on the continuum limit for both cases and derive specific closed-form solutions for the reproductive number for a configuration model static network, and an activity model temporal network.

Fully-mixed population Consider a fully-mixed population of N individuals and an infectious rate function, β(τ ). In our formalism, we denote the fraction of the population in the susceptible, jth infectious stage, and the recovered stage as S, I j , j = 1 . . . n, and R respectively and note that S + n j=1 I j + R = 1 by conservation. Assuming that an individual's infection status is independent of the infection status of its neighbors, as done in Ref. [16] , we can write the following system of mean-field equations as

By construction, an infected node will always transition through all the infectious states until it reaches the recovered state. However, we are not interested in whether infected nodes transition through all the states, but rather whether susceptible nodes become infected. In Ref. [22] , the authors introduce the notion of a next generation matrix (NGM) which decomposes the linearized system into infectious transmissions, T , and noninfectious transitions, Σ, where transmissions move susceptible nodes to infected compartments and transitions move infected nodes to other infectious states. As done in Ref. [22] , we exclude the susceptible and recovered states. The linearized system can be written as

where I = (I 1 , . . . , I n ) T . We split the matrix into transmissions and transitions and according to Ref. [22] , the reproductive number R 0 is given by ρ(−T Σ −1 ) which for the fully mixed case evaluates to

which matches the value found in Ref. [16] . This result indicates that any infectious rate function that has the same total infectiousness or exposure yields the same reproductive number, regardless of the particular function. This, however, does not hold for the time scale on which the epidemic spreads as we will see later.

Discrete category-mixed population

Now we consider a population with N individuals each of which belong to a category c i , i = 1 . . . n c . These mixing categories can encode many different characteristics such as degree-based mixing [23] , age-mixing [24] , spatial meta-population mixing [5] , mixing due to travel and many other types of mixing.

We denote the probability that sub-populations c i and c j interact with each other as p(c i , c j ) and the probability that a node belongs to category i as p(c i ). We discretize the infectious states not only by the progression of the infection, but by the category to which that individual belongs as well. This model has (n + 2)n c states: n c susceptible states, S c1 , . . . , S cn c ; nn c susceptible states, 

The linearized ODE is the following block-matrix system of equations: Splitting the matrix into transmissions and transitions, the next-generation matrix is

(4) Then, the reproductive number evaluates to

which indicates that the epidemic threshold depends both on the infectious exposure and the matrix of mixing probabilities and that these two quantities are independent.

For each case described prior, it is natural to want to take the limit as the number of infectious compartments approaches infinity and ∆τ → 0. For the fully-mixed case, the reproductive number becomes

and similarly, for category-based mixing, it is

Alternatively, we can treat τ as a continuous quantity and track the infectiousness, I(t, τ ), as a function of the overall time and how long an individual has been infected. When τ is continuous, ∆τ → 0 and the finite difference (I j−1 − I j )/∆τ in Eqns. (1c) and (3c) becomes a derivative with respect to τ . With these assumptions, our ODE model can be expressed as the transport equation with boundary conditions handling the infection and recovery. For the fully-mixed case, this is

The transport equation admits traveling wave solutions and this perspective lends physical interpretation to our model; an infected individual is transported through the infectious stages and the boundaries merely introduce new individuals into this transport process and remove recovered individuals at the other boundary. We can see this behavior in Fig. 4 for both static and temporal networks. Because our approach approximates the infectious rate function with discrete infectious compartments, we perform numerical experiments to analyze the number of states at which we can expect the mean-field ODE model to reasonably approximate the continuous rate function. For a small number of states, the discretized values of the infectious rate function fluctuate, leading to nonmonotone and non-smooth trends, so we only look at the viral load model with greater than 4 infectious states. As the number of infectious states is increased, the epidemic dynamics converge to that of the continuous VL model with a continuous infectious rate function. From Fig. 2 , approximately 100 infectious states are necessary to capture key features of the epidemic response.

In the following, we apply our category-mixing framework to two cases, a static degree-based configuration model and a temporal activity-based model.

Consider a network of size N with a degree sequence k = (k 1 , . . . , k N ) T and nodes connected by links at random, which specifies the configuration model (described more in Ref. [25] ). Networks generated with the configuration model may have a non-negligible number of self-loops and multi-edges in the infinite size limit [26] , leading to correlated simple networks. In this study, however, we consider a bounded degree distribution and so we can assume the configuration model to be uncorrelated for large enough N . For the standard SIR model on a configuration model network, the reproductive num-ber is R 0 = β k 2 /(γ k ) [1] . We assume that a node's degree completely specifies its dynamic behavior, which ignores effects from a node's other characteristics. From the degree sequence k, we can compute the discrete probability distribution p(k) = N (k)/N , where N (k) is the number of nodes in the degree sequence that have degree k, and the list of unique degrees in the degree sequence, k u . From our general formalism in Section III B, the degree mixing matrix is

where k u p = (k 1 p(1), . . . , k max p(k max )) T and k u = (k 1 , . . . , k max ) T . The largest eigenvalue of this matrix is k 2 / k and so the reproductive number is

Setting γ = 1/τ R and β = β(τ ) = τ R 0 β(τ )dτ /τ R for the SIR model yields the reproductive numbers derived in Ref. [1] .

Our category-based framework applies not only to static contact structures, but to temporal networks as well. We consider the activity model first presented in Ref. [27] . Given a temporal network of size N , suppose that each node i has an activity rate a i , which denotes the probability per unit time that the node is active. At each discrete time, each node is either active or idle, and each active node forms m connections with other nodes, active or inactive. Unlike degrees which are discrete for an unweighted network, these activity rates are continuous, and to use our category-based mixing framework, we assume that we can bin these rates into discrete categories, a i , i = 1 . . . n a and later take the continuum limit as before. We denote the probability that a node has an activity rate a i as p(a i ). Then the probability that nodes with activity rates a i and a j are connected at any given time is (a i + a j ) m N and the time-averaged mixing matrix is

which can be written P = 1b T + cp T where b = (m a 1 p(a 1 ), . . . , m a na p(a na )) T , c = (m a 1 , . . . , ma na ) T , and p = (p(a 1 ), . . . , p(a na )) T . Observing that this is a rank-2 matrix, the analytical solution for the Perron-Frobenius eigenvalue is (m a + m a 2 ) and

In Ref. [27] , they derive the epidemic threshold for the activity model as β/γ = 2 a /( a + a 2 ). As before, setting γ = 1/τ R and β = k β(τ ) = 2m a β(τ ) yields the same result.

In Ref. [28] , the authors consider heterogeneous susceptibility and recovery rate for the SIR model. Similarly, we now relax the assumption that the infectious rate function is the same for every individual. We extend our results in Section III B for a distribution of infectious rate functions over the population. In our analysis, we assume that the particular infectious rate function is distributed independently of any other nodal characteristic such as its degree. We denote p b (b) as the fraction of the population with an infectious rate function of β b (τ ) and an associated recovery time of τ R b , where the number of unique infectious rate functions is n b . We enforce that the number of infectious states regardless of recovery time is n so the time between infectious compartments is n∆τ b = τ R b . We define the discretized values β i (τ j ) = β i (j∆τ i ) as β i j and denote the jth infectious stage with infectious rate function β b (τ ) and category c as I b,c j . Then the mean-field equations become

Linearizing these equations, we obtain I = AI, where A = Σ + T . Σ and T are each n × n block matrices of size n c n b × n c n b with blocks of size n c × n c .

Then, the reproductive number (the maximal eigenvalue of −T Σ −1 ) is

As n → ∞, ∆τ b → 0 for every b and we obtain

which is the value obtained for the category-mixed case with the key difference that the exposure is now the average exposure with respect to the distribution of infectious rate functions.

We compare the time dynamics of the SIR model with that of the VL model with different infectious rate functions. For the following figures, we fixed N = 10 4 , R 0 = 3, τ R = 21 days, and arg max τ β(τ ) = 4 days, unless otherwise noted. We considered the configuration model with a power-law degree distribution p(k) ∝ k −3 on [10, 1000] and the activity model with activity rates p(a) ∝ k −3 on [0.01, 1], m = 10, and ∆t = 1. We used the following contagion models: the VL model with β Γ (τ ) ∝ τ exp(τ /4) as in Ref. [12] , the VL model with a constantvalued infectious rate function, β const (τ ) = β Γ (τ ) , and the SIR model with a single infectious rate of β = β Γ (τ ) for the configuration model and β = 2m a β Γ (τ ) for the activity model. These relations were chosen such that the reproductive numbers are the same for each infection model.

We simulated all the contagion models described in discrete time with ∆t = ∆τ = 1. We simulated the SIR model as a discrete time Markov process using the parameters γ = 1/τ R and β = β Γ (τ ) and β = 2m a β Γ (τ ) for the configuration and activity models respectively. For the viral load model, we store the time at which node i has been infected as t * i and at time t, the rate of infection of that node is β(t − t * i ), for example, and when t−t * i ≥ τ R , the node recovers. When simulating on temporal networks, we store the temporal network as an array, where each entry is a static network corresponding to a particular snapshot in time.

From Fig. 3 we see that the peak of the SIR model is delayed relative to both viral load models and for the static case, the epidemic peak is significantly less pronounced. We comment that the viral load model fundamentally changes the time scale of the epidemic when compared to the SIR model. Not all infectious compartments are created equal, however; someone at their peak infectiousness contributes much more to the spread of an epidemic than someone who has just gotten infected or almost recovered. For this reason, we now plot the number of individuals in each infectious stage over time. We now relax the assumption that β Γ (τ ) and τ R are identical for each member of the population. We assume that arg max τ β(τ ) ∼ Uniform(2, 6) and that τ R ∼ Uniform(16, 26) similar to Ref. [10] . At given times t, we plot the number of individuals as a function of the infectious duration τ i = t − t * i and t. We see traveling wave behavior as described in Section III C for both static and temporal networks. The amplitude of this wave varies in response to the introduction of new infected individuals, but the distribution shows a clear transition to the latter infectious stages as the epidemic progresses. This behavior is corroborated by the three normalized vertical cross-sections, showing the probability distribution at selected times. We notice that, despite identical values of β i (τ ) and τ Ri for every node, the temporal behavior is different for static and temporal networks. For the temporal network, it seems evident that individuals with the longest infection duration seem to be driving the epidemic based on the minimal decrease in individuals for large τ in comparison to the static network case.

We also plot the epidemic extent for different values of R 0 in Fig. 5 to validate our predictions of the reproductive number. For each data point, we averaged over 100 simulations but use the same network realization for all simulations for both the configuration and activity models. We ran each simulation until there were no longer any infected individuals.

We see that for both static and temporal networks, the predictions from our theory do as well as the predictions for the SIR model in Refs. [1] and [27] . The gradual transition is due to the heterogeneity of the networks and agrees with prior results on power-law networks [28] .

In our analysis, we theoretically derived and numerically validated predictions of the population reproductive number for static and temporal networks for a contagion model accounting for time-dependent infectiousness. We see that time-dependent infectiousness causes a funda-mental change in the time dynamics compared to the dynamics of the SIR model, despite an epidemic threshold matching classical theory.

Although time-dependent infectiousness does not affect predictions on whether an epidemic will initially grow or die out, it has strong implications how the epidemic progresses in time. In the continuum limit, the viral load model can be written as the transport equation PDE with an infectious boundary condition, which indicates that distribution of τ progresses in time like a traveling wave and this prediction is validated with numerical simulations.

In this study, we have only considered the population reproductive number, though it is well known that merely studying the population reproductive number without examining the heterogeneity in the number of secondary infections leaves out key information [8] . Superspreading events are the result of this stochasticity and can often be responsible for the transmission of an epidemic. The VL framework could be used to model the distribution of secondary infections resulting from a combination of contact-based and infectiousness-based heterogeneity.

Epidemic spreading in correlated complex networks

Epidemic dynamics in finite size scale-free networks

Spread of epidemic disease on networks

Epidemic processes in complex networks

Modeling the spatial spread of infectious diseases: The GLobal Epidemic and Mobility computational model

Modeling the Spatiotemporal Epidemic Spreading of COVID-19 and the Impact of Mobility and Social Distancing Interventions

A model for the spread of an epidemic from local to global: A case study of COVID-19 in India

Superspreading events in the transmission dynamics of SARS-CoV-2: Opportunities for interventions and control

Social confinement and mesoscopic localization of epidemics on networks

Test sensitivity is secondary to frequency and turnaround time for COVID-19 surveillance, medRxiv

Transmission of COVID-19 in 282 clusters in Catalonia, Spain: A cohort study

Temporal dynamics in viral shedding and transmissibility of COVID-19

Realistic Distributions of Infectious Periods in Epidemic Models: Changing Patterns of Persistence and Dynamics

Mathematics of Epidemics on Networks: From Exact to Approximate Models

Generality of the Final Size Formula for an Epidemic of a Newly Invading Infectious Disease

The differential infectivity and staged progression models for the transmission of HIV

Message passing approach for general epidemic models

Mean-field models for non-Markovian epidemics on networks

Pairwise approximation for SIR-type network epidemics with non-Markovian recovery

Rethinking Covid-19 Test Sensitivity -A Strategy for Containment

Dynamics of Multi-stage Infections on Networks

The construction of next-generation matrices for compartmental epidemic models

Edge-based compartmental modelling for infectious disease spread

Inferring highresolution human mixing patterns for disease modeling

Configuring Random Graph Models with Fixed Degree Sequences

Generation of uncorrelated random scale-free networks

Activity driven modeling of time varying networks

How heterogeneous susceptibility and recovery rates affect the spread of epidemics on networks

I would like to thank Ren Stengel for working on stochastic simulations that were helpful in framing this study, Daniel Larremore for many helpful conversations and theoretical insights, Juan G. Restrepo for helpful suggestions and draft edits, and Subekshya Bidari for her help in formulating the PDE model.

All code used in this study can be found at https://github.com/nwlandry/time-dependentinfectiousness [29] .