Who is the infector? General multi-type epidemics and real-time susceptibility processes
Tom Britton Ka Yin Leung Pieter Trapman
math.PR 92D30, 60K35



  
   
                
   
  
theoremTheorem[section]
corollary[theorem]Corollary
lemma[theorem]Lemma
proposition[theorem]Proposition
remark[theorem]Remark
conjecture[theorem]Conjecture
assumption[theorem]Assumption













Who is the infector? General multi-type epidemics and real-time susceptibility processes
    Tom Britton^1, Ka Yin Leung^1 and Pieter Trapman^1
    December 30, 2023
========================================================================================




We couple a  multi-type stochastic epidemic process with a directed random graph, where edges have random lengths.
This random graph representation is used to characterise the fractions of individuals infected by the different types of vertices among all infected individuals in the large population limit.   
For this characterisation we rely on theory of multi-type real-time branching processes.
We identify a special case of the two-type model, in which the fraction of individuals of a certain type infected by individuals of the same type, is maximised among all two-type epidemics approximated by branching processes with the same mean offspring matrix. 


Keywords: Epidemics, Multi-type branching process approximation, Susceptibility processes, Directed random graphs


2010 Math. Subj. Classification:
Primary 92D30; 60K35; 
Secondary 05C80; 60J80

[1]Department of Mathematics, Stockholm University, 106 91 Stockholm, Sweden.
 Email: {tom.britton, kayin.leung, ptrapman}@math.su.se


§ INTRODUCTION

Mathematical models have proven to be successful in understanding infectious disease dynamics. Often, the focus is on (controlling) the beginning of an epidemic. In many of those models branching process approximations and the concept of the basic reproduction number R_0 (which corresponds to the offspring mean of the approximating branching process) play an important role. In the current paper we consider  the entire epidemic outbreak instead. We do so in a setting where there are multiple types of infected individuals. Suppose that a large epidemic outbreak has taken place in the population. Then a certain fraction of the population will have become infected. The question that we concern ourselves with is: what fraction of the infected population of a certain type j was infected by individuals of type i? Here i,j=1,… K, where K is the number of types in the population. In other words, who is the infector?

This question is not so straightforward to answer. Timing of events plays an essential role. We approach this question using an epidemic graph construction, which is used as a tool in proving the two main theorems of this paper. But, the construction proves to be interesting in itself, and a substantial part of this paper is devoted to this construction. Using the graph representation of the epidemic, we consider susceptibility processes. Susceptibility sets were introduced in infectious disease modelling by Ball and co-authors <cit.>. The susceptibility set of a vertex v consists of all vertices u in the vertex set from which there is a path from u to v in the (restricted) epidemic random graph. The epidemic process and epidemic random graph can be coupled in such a way that u is in the susceptibility set of v if and only if v is infected during the epidemic conditioned on u being initially infectious. Susceptibility sets have proven to be important tools in proving results concerning e.g. the final size of different epidemic models <cit.> (see also <cit.>). However, for the research question in this paper we need to consider the susceptibility process instead, i.e. we need to take timing into account. 

In this paper, we also consider a multi-type (backward) branching process. This branching process is constructed in such a way that the distribution of the tree-like graph corresponding to it, is the same as the susceptibility process (up to a certain time). The coupling between the susceptibility process and the branching process enables us to prove the main results of this paper formulated in Theorems <ref> and <ref>. Using existing theory for branching processes we find an answer in Theorem <ref> for the question “who is the infector?” by means of an expression for the expected fraction ρ_ij of infected individuals of type j that are infected by individuals of type i, conditioned on that there is a large outbreak in a population (as the population size tends to infinity), i,j=1,… K. 
In general, this expression remains rather implicit. However, for a special class of models, we are able to obtain upper and lower bounds for the quantities ρ_ij of interest (Theorem <ref>) if we keep the fractions of individuals of the different types and the expected number of infectious contacts between different types of individuals fixed. As the ρ_ij of Theorem <ref> are rather implicit, these bounds allow us to gain more insights in the importance of different types of infected individuals in the transmission dynamics in the population. 

The class of models that allows for identifying upper and lower bounds are the topic of interest of our twin paper <cit.>. This paper <cit.> is motivated by infectious diseases, such as influenza and chlamydia, for which we can categorise infected individuals as symptomatic or asymptomatic (showing no apparent signs of the disease) giving rise to two types of infected individuals. Asymptomatically infected individuals are generally hard to detect by public health authorities. Therefore, we would like to gain insights in their role in the transmission process and determine whether asymptomatically infected individuals often play the role of the infector. 

The structure of this paper is as follows. In Section <ref> we introduce the model, the notation and the two main theorems of the paper. Next, in Section <ref>, we discuss the construction that enables us to prove the desired results. The proofs are then presented in Section <ref>. We end with a short discussion in Section <ref>.




§ MODEL, NOTATION AND MAIN RESULTS





 §.§ Model and notation

We denote the number of elements in a set 𝒜 by |𝒜|. For k ∈ℕ we use the notation [k] = {1,2, ⋯, k}. We use ℕ for the strictly positive integers and ℕ_0 = ℕ∪{0} for the non-negative integers. 

Unless specified otherwise, limits are for population size n →∞. We say that an event happens with high probability (w.h.p.) if the probability of the event converges to 1 as n →∞. We adhere to the usual order notation, i.e. f =O(g) means that lim sup_n →∞ |f(n)/g(n)| < ∞ and f =o(g) means that lim_n →∞ |f(n)/g(n)| =0. In addition, for a sequence of random variables {X^(n); n∈ℕ}, we write X^(n) = O_p(g) if |X^(n)/g(n)| is bounded in probability and X^(n)=o_p(g) if X^(n)/g(n) → 0 in probability. See <cit.> for a discussion of this notation.  

We consider a population of n individuals where V^(n) denotes the set of all individuals. For some of our results we consider a sequence of models in growing populations, i.e. for n →∞. If no confusion is possible we write V= V^(n). We assume that there are K types of individuals. For i ∈ [K], let V_i = V_i^(n) be the set of vertices of type i and n_i = |V_i| the number of vertices of type i. 




Within this sequence of populations we consider the spread of an infectious disease in which individuals are either susceptible or infected. When individual v ∈ V_i^(n) is infected it makes contacts to different individuals of type j (i,j ∈ [K]) according to a point process ξ_v^j = {ξ_v^j(t);t≥ 0} on the interval [0,∞). If there are more than n_j points in the process, then only the first n_j points represent contacts. The time parameter in the definition of ξ_v^j represents the time since infection of individual v. The processes ξ_v^1, ξ_v^2, ⋯ and  ξ_v^K may be dependent, and their joint distribution may depend on the type i of v. However, the point processes associated to different individuals are independent, i.e. the vectors of processes {(ξ_v^1, ⋯ξ_v^K);v ∈ V} are independent. For convenience we introduce the vector of stochastic process (ξ_i1, ξ_i2,⋯, ξ_iK), which is distributed as (ξ_v^1, ⋯ξ_v^K) for v ∈ V_i, i ∈ [K]. We make the following assumptions.


n^-1 n_i → p_i>0  for all i ∈ [K]. 
Furthermore, 

    max_i,j∈ [K]|p_i/p_j-n_i/n_j|=O(1/n).




For all i ∈ [K], the distribution of (ξ_i1, ξ_i2,⋯, ξ_iK) is independent of population size n and  for all i,j ∈ [K],

    m_ij=𝔼[ξ_ij(∞)]<∞,

i.e. the expected number m_ij of contacts that an individual of type i makes with individuals of type j  is finite, for all i,j ∈ [K].  Furthermore, we assume that |(ξ_i1, ξ_i2,⋯, ξ_iK)| does a.s. only contain jumps of size 1 and its distribution has no atoms. 



There exists a constant κ <1 such that for all i ∈ [K], 
    ℙ(max_v ∈ Vξ_v^i(∞) < n^κ) → 1          n →∞.



For future reference, we let M={m_ij}_i,j∈[K] denote the matrix with the m_ij (as defined in Assumption <ref>) as elements.
Assumption <ref> guarantees that w.h.p. all non-trivial paths in the epidemic graph defined below have different lengths. 
Note that Assumption <ref> is easily met, e.g. if for all v ∈ V and i ∈ [K] there exists ϵ>0 and ℓ_0 ∈ (0,∞), for which ℙ[ξ_v^i(∞)>ℓ] < ℓ^-(1+ϵ) for all ℓ>ℓ_0, then for κ∈ (1/(1+ϵ),1)

    ℙ(max_v ∈ Vξ_v^i(∞) < n^κ)    =  ∏_v ∈ V [1-ℙ( ξ_v^i(∞) > n^κ)] 
       ≥ 1 - ∑_v ∈ Vℙ( ξ_v^i(∞) > n^κ) 
       ≥ 1 - n max_j ∈ [K]ℙ( ξ_ji(∞) > n^κ)
        >1 -n × n^-(1+ϵ)κ
       → 1.


At the points of ξ_v^j (j ∈ [K]), v contacts an individual from V_j. The individuals that are contacted are uniformly chosen without replacement. If v is of type j, we allow for v to be among the contacted individuals in ξ_v^j. If the individual that is contacted is still susceptible at the time of the contact then it becomes infected. 

Note that we may assign the point processes {ξ_v^j;v ∈ V, j ∈ [K]}  and decorate the points with the labels of the individuals these points represent contacts to, already before the epidemic starts. This allows us to create a new set of random variables {η(u,v);u,v ∈ V}. Assume that individual v is of type j ∈ [K]. If there is a point in ξ_u^j with label v, then η(u,v) takes the value of this point. If there is no such point in ξ_u^j, set η(u,v)= ∞. Observe that the distribution of {η(u,v);u,v ∈ V} depends on n in this construction, because the probability that v is chosen as a label is decreasing in n. Note that there is a broad class of models that satisfy  Assumptions <ref>-<ref>, see Remark <ref> for an example of a specific class of epidemic models.

This construction provides us with a graph representation of the population and the epidemic on it. We construct the weighted random graph G=(V,E) as follows. The edge set E consists of all directed pairs (u,v) ∈ V × V, with u ≠ v. For edge (u,v) ∈ V × V, we say that u is the tail of (u,v) and v is its head. The weight of edge (u,v) is given by η(u,v) for all (u,v) ∈ E. For some of our arguments we restrict to the weighted edge set E' ⊂ E of all edges (u,v) ∈ E with finite weight, i.e. (u,v) ∈ E' if and only if  η(u,v)< ∞. The corresponding random graph is denoted by G'=(V, E'). For i,j ∈ [K],  u ∈ V_i and v ∈ V_j, we also introduce the random variable η_ij, which is distributed as η(u,v) conditioned on η(u,v)<∞ (i.e. conditioned on (u,v) ∈ E').
Observe that 

    ℙ(η_ij≤ t)    = ℙ(η(u, v) ≤ t|η(u, v) < ∞) 
       = ℙ(η(u,v) ≤ t)/ℙ(η(u,v) < ∞) = 1/n_j𝔼[ξ_ij(t)]/1/n_j𝔼[ξ_ij(∞)]= 𝔼[ξ_ij(t)]/m_ij,

which is independent of n. We mainly consider the random graph G'=(V,E'). See Figure <ref> for an example of G'.








An ordered set of distinct vertices π =(v_1, v_2, ⋯, v_m) is a path in G' if (v_i,v_i+1) ∈ E' for all integers i ∈ [m-1]. With some abuse of terminology, we sometimes refer to the set of edges connecting the vertices of π as the path π and speak of a path in E'. The length of a path π =(v_1, v_2, ⋯, v_m) is ℓ(π) = ∑_i=1^m-1η(v_i,v_i+1). Assumption <ref> implies that all non-zero paths in E' have different lengths with probability 1. Let Π_uv be the set of all paths from u to v in E' and define the (quasi) distance from u to v as d(u,v)=min_π∈Π_uvℓ(π). As an example, in Figure <ref> the distance from a to d is given by 
    d(a,d) =min(η(a,b)+η(b,d), η(a,c)+η(c,d))=1.8.

In general, d(u,v)≠ d(v,u) if u≠ v since the graph G' is directed.   Therefore d is actually a quasi-distance. We say that d(u,v)=∞ if Π_uv = ∅. Furthermore, d(v,v)=0 for all v∈ V. 

Next we formulate an “irreducibility” assumption:


For every i,j∈[K], there is w.h.p. a path from a vertex in V_i to a vertex in V_j in E'.
 








An epidemic process is reproduced from G' as follows. 
Let V_init=V_init^(n) be the set of vertices that are initially infected. This set may be predetermined or randomly selected and satisfies the following assumption.


The set V_init^(n) of vertices that are initially infected satisfies |V_init^(n)|=O_p(1).
 

We set σ_v to be the time between the start of an epidemic until v becomes infected, v ∈ V, i.e. σ_v = inf{d(u,v);u ∈ V_init}. In particular, σ_v =0 for v ∈ V_init.
So, suppose we let V_init={a} in Figure <ref>, then σ_a=0, σ_b=0.3, σ_c = 1.3,  σ_d = 1.8 etc.
Moreover, if σ_v= ∞, then v will never become infected. Note that G' contains some redundant information regarding the epidemic: (i) there are edges with finite weights with infection time of the tail being ∞, i.e. the tails of those edges will never get infected, (ii) if σ_v < σ_u + η(u,v), then the edge (u,v) represents a contact between two already infected individuals. 


The framework of this section includes the SEIR (Susceptible → Exposed → Infectious → Recovered) epidemics as follows. In the SEIR framework, exposed, infectious and recovered individuals are all counted as infected. Assign to every vertex v ∈ V a random latent period L_v and a random infectious period ι_v, which might be dependent on L_v. The vectors {(L_v,ι_v);v ∈ V} are independent and their distribution functions only depend on the type of the vertex. Suppose that v is of type i. Conditioned on (L_v,ι_v), let the processes {ξ̂_v^j; j ∈ [K]} be independent homogeneous Poisson processes on the interval (L_v, L_v + ι_v) with intensity p_j λ_ij. At the points of this point process, v makes contacts to  vertices in V_j chosen uniformly with replacement. By keeping only the points in {ξ̂_v^j; j ∈ [K]} whose label did not appear before in this process, we obtain {ξ_v^j; j ∈ [K]}. If L_v and ι_v  are exponentially distributed, then the SEIR epidemic is a Markov process that is often referred to as the Markov SEIR epidemic.

At time t, a vertex v is susceptible if t< σ_v, Exposed if t ∈ [σ_v, σ_v + L_v), while it is Infectious if t ∈ [σ_v + L_v, σ_v + L_v+ ι_v). Finally v is Recovered if t ≥σ_v + L_v+ ι_v.

If ℙ(L_v =0)=1, then we are in the so-called SIR (Susceptible → Infectious → Recovered) epidemics framework. In addition, if ι_v is exponentially distributed then the process is called a Markov SIR epidemic. 


Finally, to conclude this section, we introduce the basic reproduction number R_0. As this is possibly the most studied quantity in mathematical modelling of the spread of infectious diseases, no epidemic modelling paper would be complete without at least mentioning R_0. In a single-type epidemic model R_0 is defined as the expected number of infectious contacts made by a newly infected individual in an otherwise susceptible population. The multi-type equivalent of R_0 is given by the dominant eigenvalue of the K × K matrix M={m_ij}_i,j∈[K] (see Assumption <ref>). R_0 is always real and strictly positive <cit.>,<cit.>. We say that the epidemic process is supercritical (resp. critical, resp. subcritical) if R_0>1 (resp. R_0=1, resp. R_0<1). As n →∞, the epidemic becomes large  with positive probability if and only if the process is supercritical <cit.>. In the remainder of the paper we assume that R_0>1.















 



 §.§ Main results

Because we assume that R_0>1, a large outbreak occurs with positive probability and at the end of such a large outbreak, a fraction of the infected individuals is of type j w.h.p. However, the question that we are interested in is: what fraction of those individuals was infected by individuals of type i, for i,j∈[K]?
So, if we denote set of infected individuals of type j by ℐ_j, and the subset of those individuals, which are infected by individuals of type i  by ℐ_j^i, then we are interested in |ℐ_j^i|/|ℐ_j|. 
For example, in Figure <ref>, if a represents the initially infectious individual,
then there are two individuals of type 1 that are infected by an individual of type i (namely the individuals represented by vertices c and g) and one individual of type 1 that is infected by an individual of type 2 (namely the individual represented by vertex d).
So, the fraction of individuals infected by an individual of type 1 among all infected  individuals of type 1 is 2/3 (we exclude vertex a because it was initially infected, and not infected by another individual in the population).
The  question posed above leads to the main results that are formulated in Theorems <ref> and <ref> below. For our results, we need to define a multi-type branching process {𝒵(t); t ≥ 0}. The branching process is defined as follows. Let Z^i(t) be the number of particles in the branching process at time t if the process is started by a single  particle of type i, i ∈ [K]. For i,j ∈ [K], particles of type j in 𝒵(t) give birth to particles of type i according to a Poisson process with intensity p_i/p_j𝔼[ξ_ij(da)], where a is the age of the particle of type j. All of these Poisson processes are independent.

The definition of the branching process 𝒵(t) is such that we can apply theory from <cit.>. In particular, there exists a Malthusian parameter α>0 and a random variable W^i, such that 

    e^-α tZ^i(t) → W^i
 
and 

    ℙ(W^i ∈ (0, ∞))=1 -ℙ(W^i =0) = ℙ(Z^i(t) →∞).
 
Let W^i(r) for r ∈ℕ be independent copies of W^i. We now state our first main result. 



Conditioned on the occurrence of a large outbreak,
the fraction of infected individuals of type j that are infected by an individual of type i during an outbreak in a population of size n converges in probability to ρ_ij. Here ∑_i=1^Kρ_ij=1, and

    ρ_ij=1/ℙ(Z^j(t) →∞)𝔼[∑_r=1^X_ije^-ατ_ijrW^i(r)/∑_k =1^K ∑_r=1^X_kj e^-ατ_kjrW^k(r)(∑_k =1^K ∑_r=1^X_kjW^k(r) >0)
    ],

where, for k ∈ [K] and r ∈ℕ, the random variables τ_kjr are independent with distribution function ℙ(τ_kjr <a) = 𝔼[ξ_kj(a)]/m_kj and X_kj is Poisson distributed with expectation p_k/p_jm_kj.



Note that, in general, it is hard to give a more explicit expression for ρ_ij than (<ref>). Often, there is no explicit description of the distribution of W^k(r). We are able to obtain bounds for (<ref>) for the important special case discussed in <cit.>, leading to Theorem <ref>. 

The model of <cit.> is as follows. We consider K=2, and (ξ_v^1,ξ_v^2), v ∈ V, obtained from a single marked point process ξ_v. In this process ξ_v the points get independently mark 1 with probability p_1 and mark 2 otherwise. Then the process ξ_v^1 consists of the points with mark 1, while ξ_v^2 consists of the points with mark 2. By construction of the process, the probability that an ultimately infected vertex of type i is infected by a  vertex of type 1 is the same for i=1 and i=2, i.e. ρ_11=ρ_12=ρ_1. With some abuse of notation we write ξ_i for a point process distributed as ξ_v for v ∈ V_i, i ∈ [2]. Furthermore, note that m_ij=𝔼[ξ_ij(∞)]=p_j𝔼[ξ_i(∞)], i.e. we can write m_ij=p_j m̃_i with m̃_i=𝔼[ξ_i(∞)]. Here m̃_i can be interpreted as the expected number of secondary cases caused by a newly infected individual of type i in an otherwise susceptible population (one can think of the m̃_i as the type-specific reproduction numbers). Then the basic reproduction number is R_0=p_1m̃_1+p_2m̃_2, with p_2=1-p_1. Indeed, a newly infected individual is of type i with probability p_i and the expected number of secondary cases it produces is m̃_i, i=1,2.  


For this model we can compute ρ_1^- and ρ_1^+, the minimum and maximum fraction of the infected vertices that are infected by type 1 vertices, if the matrix M and p_i are held fixed for i ∈ [K]. We let q̃_1 be the smallest positive solution in (0,1] of 

    x= e^-p_1 m̃_1(1-x)
 
and q̃_2 the smallest positive solution in (0,1] of

    x= e^-p_2 m̃_2(1-x).
 
Furthermore, we let q be the unique solution in (0,1) of 

    x  =  e^-(1-x)(p_1 m̃_1+ p_2 m̃_2)=e^-(1-x)R_0.

Note that we assume that R_0=p_1m̃_1+ p_2m̃_2>1, so the unique solution q∈(0,1) exists. We can interpret 1-q̃_1 (resp. 1-q̃_2) as the final fraction of the population that ultimately gets infected when only individuals of type 1 (resp. type 2) are able to transmit, conditional on a large outbreak. Furthermore, 1-q can be interpreted as the final fraction of the population that ultimately gets infected, conditional on a large outbreak (or conversely, q is the fraction of the population that remains susceptible throughout the epidemic). 



Consider the two-type model described above. In the limit as population size n→∞ and for fixed p_1, m̃_1 and m̃_2, the fraction of ultimately infected vertices that is infected by type 1 vertices is bounded from above by 

    ρ_1^+ = ( 1 -  p_1 m̃_1 (q̃_1 + q)/2) q̃_1-q/1-q.



By interchanging the role of the types 1 and 2, we also obtain the lower bound ρ_1^-. Indeed, note that ρ_1^-=1-ρ_2^+, where

    ρ_2^+ = ( 1 - p_2 m̃_2 (q̃_2 + q)/2) q̃_2-q/1-q.

In other words, for any  point process {(ξ_v^1,ξ_v^2); v∈ V} that satisfies the assumptions of Section <ref> and that can be obtained from independently marking points of a one-dimensional point process, the fraction ρ_1 of infected individuals that are infected by individuals of type 1 is bounded by ρ_1^- and ρ_1^+, i.e. ρ_1^-≤ρ_1 ≤ρ_1^+. 


In Section <ref> we consider a more general setting than the one in Theorem <ref>. Instead of assuming a single marked process ξ_i, one can consider general distributions (ξ_i1, ξ_i2) and obtain bounds (<ref>) for ρ_21^- and (<ref>) for ρ_11^+ (and, by interchanging the roles of types 1 and 2, bounds ρ_12^- and ρ_22^+). As this is somewhat more involved, we choose to present the bounds in the form of Theorem <ref> instead.




§ SUSCEPTIBILITY PROCESS, BACKWARD BRANCHING PROCESS AND THE COUPLING


Throughout we assume that all random variables and processes are defined on a suitable rich enough probability space, which we do not specify.



 §.§ The susceptibility process

In this subsection we use the idea of susceptibility sets <cit.>, and construct this set through a stochastic process: the susceptibility process.

We define the susceptibility process {𝒮^(n)_v(t);t ≥ 0} as 

    𝒮^(n)_v(t) = {u ∈ V; d(u,v) ≤  t}.
  
Note that 𝒮^(n)_v(t) is non-decreasing in t. The susceptibility set of vertex v is defined as 𝒮^(n)_v = lim_t →∞𝒮^(n)_v(t), i.e. the susceptibility set S_v^(n) of v consists of all vertices u∈ V for which there is a path from u to v in G'. As an illustration, in Figure <ref>, the susceptibility set of vertex f  is given by  𝒮^(n)_f = {a,b,c,d}, while 𝒮^(n)_f(t=1.1) = {c,d}. Note that 𝒮^(n)_v ∩ V_ init = ∅ if and only if v remains uninfected throughout the epidemic, i.e. if and only if there is no path in G' from V_ init to v. Also note that the susceptibility set may contain vertices of different types. 

Define for a ≥ 0 and j ∈ [K], 

    𝒮^(n)_v(t;a,j) = {u ∈𝒮^(n)_v(t)∩ V_j; d(u,v) > t- a}.

That is, 𝒮_v(t;a,j)  consists of the  vertices of  type j in 𝒮^(n)_v(t), that are not yet part of 𝒮^(n)_v(t-a).  


Let v ∈ V ∖ V_init be a randomly chosen vertex. We derive the susceptibility process {𝒮^(n)_v(t);t ∈ (0,t_*)} by constructing part of the random graph G' around vertex v by means of an exploration process {Ĝ(ℓ)}= {Ĝ(ℓ); ℓ∈ℕ_0} in which vertices in the susceptibility process are explored one at a time. We note that {Ĝ(ℓ)} depends on v. Here t^*= t^*(n) is a given time (we defer the specification of t^* until (<ref>) below). The process {Ĝ(ℓ)} allows us to couple {𝒮^(n)_v(t)} with an appropriate branching process. In this way, we can make the coupling between the susceptibility process and the backward branching process that is needed to prove Theorems <ref> and <ref> in Section <ref>.

Before we define the exploration process {Ĝ(ℓ)} around vertex v, we introduce some additional variables and terminology. Ĝ(ℓ) is a 4-tuple: 

    {Ĝ(ℓ)= (V̂^a(ℓ),V̂^p(ℓ),V̂^e(ℓ),Ê(ℓ));ℓ∈ℕ_0}.

Here Ê(ℓ) denotes the edge set of Ĝ(ℓ). Vertices in Ĝ(ℓ) can be `active', `passive', or `explored'. The sets of these vertices are denoted by V̂^a(ℓ), V̂^p(ℓ), and V̂^e(ℓ), respectively. The sets V̂^a(ℓ), V̂^p(ℓ), V̂^e(ℓ), and Ê(ℓ) are defined in the construction below (where also their names will become apparent). Finally, before we explain the construction, we mention that throughout the process we may `flag' the process (see step 4). This flagging plays a role when coupling the exploration process with a branching process to represent the susceptibility process. The construction is as follows. 


  1) 
To set the initial conditions of the construction, let V̂^a(0) = v, V̂^e(0) = ∅, and Ê(0) be the set of all edges in E' for which v is the tail, whereas V̂^p(0) is the set of all heads of edges in Ê(0) that are in V∖ v.


  2) 
For ℓ∈ℕ_0, assume that V̂^a(ℓ)≠∅ and that there exists a vertex in V̂^a(ℓ) from which there is a path in Ĝ(ℓ) to v of length at most t^*. In step ℓ+1 pick (according to some rule) one of the vertices from V̂^a(ℓ), from which there is a path in Ĝ(ℓ) to v of length at most t^*. Say that this vertex is v' ∈ V_j. Move v' to the set of explored vertices, i.e. V̂^e(ℓ+1) = V̂^e(ℓ) ∪ v'. Assign to v' independently a binomial random number x(v';i) with parameters n_i and m_ij/n_j.


  3) 

The remainder of step ℓ+1 is split up  in ∑_i=1^K x(v';i) sub-steps as follows. We introduce

    {G^*(ℓ,ℓ');ℓ∈ℕ_0, ℓ' ∈{0,1,⋯, ∑_i=1^K x(v';i)}},

where

    G^*(ℓ,ℓ')= (V^*,a(ℓ,ℓ'),V̂^*,p(ℓ,ℓ'),V̂^*,e(ℓ,ℓ'),Ê^*(ℓ,ℓ')).

Furthermore, set G^*(ℓ,0)=Ĝ(ℓ). Next, let ℓ' ∈(∑_i'=1^i-1 x(v';i'),∑_i'=1^i x(v';i')], where the empty sum is 0.

  4)
In the ℓ'-th sub-step of step ℓ+1 pick uniformly a vertex from V_i.
If we pick a vertex we have picked in one of the ℓ'-1 sub-steps before we say that the exploration process is flagged. 
In that case we choose new vertices from V_i until we obtain a vertex that was not chosen in the previous ℓ'-1 sub-steps. This step is equivalent to picking the vertices without replacement.
Say that the vertex that is picked is v”. 


  4a) If v”∈V̂^a(ℓ) ∪V̂^e(ℓ), then nothing changes in the exploration graph, i.e. G^*(ℓ+1,ℓ')= G^*(ℓ+1,ℓ'-1). This is because if v”∈V̂^a(ℓ) ∪V̂^e(ℓ), then we already have explored the edges with tail v”.  


  4b) If v”∈ V ∖(V̂^a(ℓ) ∪V̂^e(ℓ)), assign to v” the vector of point processes (ξ_v”^j', j' ∈ [K]). The distribution of (ξ_v”^j', j' ∈ [K]) is equal to the distribution of (ξ_ij', j' ∈ [K]), given that ξ_ij contains a vertex with label v'. Assign label v' to a uniformly chosen point in ξ_v”^j, and assign uniform labels without replacement from V_j ∖ v' to the other points in ξ_v”^j and uniform labels from V_j' to the points in ξ_v”^j' for j' ∈ [K]∖ j.

If none of the newly assigned labels correspond to vertices in V̂^e(ℓ) then E^*(ℓ+1,ℓ') contains all edges in E^*(ℓ+1,ℓ'-1) plus the edges with tail v” and heads  corresponding to the labels of the points in (ξ_v”^j', j' ∈ [K]), with the obvious edge lengths. In addition, all heads of those edges which were not in V^*,a(ℓ+1,ℓ'-1) move to V^*,p(ℓ+1,ℓ') (if they were not in that set already). Furthermore, v”∈V^*,a(ℓ+1,ℓ').

If any of the newly assigned labels correspond to vertices in V̂^e(ℓ) then we flag the process and we return to the start of step 4). This last part of step 4b) is equivalent to conditioning on the event that there are no edges with tail v” and an already explored vertex as head.


  4c) Set Ĝ(ℓ+1)=Ĝ(ℓ+1, ∑_i=1^Kx(v',i)).



  5) Continue this process by increasing ℓ until there are no active vertices having a path of length less than t^* towards v in Ĝ(ℓ). Say that ℓ^* is the smallest ℓ for which this is the case.


We note the following:


  * The edge set Ê(ℓ) is a subset of E' and contains all edges in Ĝ(ℓ). 

  * All edges in E' with tails in V̂^a(ℓ) are in Ê(ℓ), but there might still be edges in E' with heads in V̂^a(ℓ)  that are not in Ê(ℓ)). 

  * Every passive vertex is the head of an edge in Ê(ℓ), but the passive vertices themselves are not explored, and none of the edges in E' with tail in V̂^p(ℓ) are in Ê(ℓ). 

  * All edges in E' with head or tail in V̂^e(ℓ) are in Ê(ℓ). 

  * The tails of edges in Ê(ℓ) are in V̂^a(ℓ) ∪V̂^e(ℓ) and their heads are in  V̂^a(ℓ) ∪V̂^p(ℓ) ∪V̂^e(ℓ).



The construction of the exploration process yields V̂^e(ℓ^*) = 𝒮^(n)_v(t^*). Furthermore, if the process {Ĝ(ℓ)} is not flagged until step ℓ^*, then the construction of V̂^e(ℓ^*) (and the distances from the vertices in this set to v) is equivalent to constructing a branching process {𝒵̃^(n)(t);t ≥ 0} up to time t^*. In this branching process, particles of type j give birth to a binomial distributed number of particles of type i, i ∈ [K], where the parameters of the binomial distributed random variable are n_i and (n_j)^-1m_ij. The number of children of the different types of particles are independent. Furthermore, the ages of the mother particles of type j at birth of a child of type i are independent and have density 𝔼[ξ_ij(da)]/m_ij.








 §.§ The (backward) branching process


We create a multi-type branching process {𝒵(t);t ≥ 0} that can be coupled to {𝒮^(n)_v(t);t ≥ 0}. The coupling is performed in Section <ref>. The branching process {𝒵(t);t ≥ 0} is constructed in such a way that the distribution of the corresponding tree-like graph is the same as that of {𝒮^(n)_v(t);t ≥ 0} up to time t^* with t^* defined by (<ref>). We leave out some of the details in the arguments. Those details can be filled in analogous to <cit.> for related models.

Without loss of generality we assume that v ∈ V_1. The multi-type (backward) branching process is as follows. The (single) ancestor is of type 1. Particles of type j give birth to particles of type i according to a Poisson process with intensity p_i/p_j𝔼[ξ_ij(da)], where a is the age of the type j particle, i,j∈[K]. All Poisson processes are independent.

The branching process {𝒵(t);t ≥ 0} is analysed using existing theory from <cit.>. First of all, the mean offspring measure of this backward branching process is defined through 

    μ^(b)_ji(dt)= p_i/p_j𝔼[ξ_ij(dt)].

Note that 

    m^(b)_ji= ∫_0^∞μ^(b)_ji(dt)= ∫_0^∞p_i/p_j𝔼[ξ_ij(dt)]=  p_i/p_jm_ij
 
is the expected number of children of type i of a particle of type j. Here m_ij is given by (<ref>) (and the corresponding matrix is M). Let M^(b)={m^(b)_ji}_j,i ∈ [K]. Straightforward matrix theory gives that M and M^(b) have the same dominant eigenvalue R_0, which by assumption is strictly larger than 1, i.e. the branching process is supercritical). 
Define 

    m̂^(b)_ji(x)= ∫_0^∞ e^-x tp_i/p_j𝔼[ξ_ij(dt)]

and let M̂^(b)(x) be the matrix with elements m̂^(b)_ji(x). Finally, let α be such that 

    M̂^(b)_ji(α)=1.

Note that K<∞, all elements of M^(b) are finite and the dominant eigenvalue R_0 of M^(b) is real and larger than 1. Therefore, α exists and is positive. 


We define the random vector Z^i(t) as Z^i(t) = (Z_1^i(t),Z_2^i(t),⋯,Z_K^i(t)), where Z_j^i(t) is the number of particles of type j ∈ [K] in 𝒵(t) if the process starts with one newborn particle of type i. With some abuse of notation, let σ(x) be the time of birth of particle x in the branching process. Note that particles in the branching process 𝒵(t) never die.  

We know that there exists an α >0 such that 

    e^-α t Z^i(t) → W^i =(W^i_1, W^i_2, ⋯ W^i_K)     ,

where W^i is a random vector that has, with probability 1, strictly positive elements on the set ∑_j=1^K Z^i_j(t) →∞ as t →∞ <cit.>.

For particle x ∈{𝒵(t);t ≥ 0}, define 

    φ_x(t) = φ_x(t;a) = (t-σ(x)< a).

Let Ẑ_j^i(t; a) = ∑_x  φ_x(t;a), where the sum is taken over all particles of type j in {𝒵(t);t ≥ 0}, i.e. Ẑ_j^i(t; a) is the number of particles of type j that have age less than a in the branching process at time t. Note that Ẑ_j^i(t; a) is increasing in a.

If Z_j^i(t) →∞ for all j ∈ [K] (i.e. if new particles keep on being born in the branching process), then

    Ẑ_j^i(t; a)/∑_j=1^K Z_j^i(t)→ c(a,j)

where c(a,j) is a constant independent of i (Theorem 2.7 of <cit.>). For our purposes we do not need to specify c(a,j) further. 

Although the branching process is not dependent on n, we want to have some bound on the number of vertices born in the branching process as a function of n. This is used in the coupling in Section <ref>. We set 

    t^*= 1-κ/4log[n]/α
 
(where κ is as in Assumption <ref>). We obtain by (<ref>) that, for all i,j ∈ [K],  

    e^-α t^*Z^i_j(t^*) = n^-(1-κ)/4 Z^i_j(t^*) → W^i_j ∈ (0,∞)

on the survival set of {𝒵(t);t ≥ 0}. In particular, this implies that

    Z^i_j(t^*)  =o(n^(1-κ)/3)





 §.§ The coupling

Note that the only difference between the branching process {𝒵̃^(n)(t)} associated to the susceptibility process of Section <ref> and the multi-type branching process {𝒵(t)} of Section <ref> is the distribution of the number of particles of type i  a particle of type j gives birth to (binomially distributed with parameters n_i and m_ij/n_j and Poisson distributed with expectation p_j/p_im_ij, resp.)





We know from <cit.> (see also <cit.>) that the total variation distance between a binomial distributed random variable with parameters n_i and 1/n_in_i/n_j m_ij and a Poisson random variable with parameter n_i/n_j  m_ij is O(1/n_i). Moreover, the total variation distance between a Poisson distributed random variable with parameter n_i/n_j  m_ij and a Poisson distributed random variable with parameter p_i/p_j m_ij is O(√(|n_i/n_j-p_i/p_j|)) =  O(1/√(n)) <cit.> (see also <cit.>). Here we have also used Assumption <ref>. By the triangle inequality this implies that the total variation distance between a Poisson distributed random variable with expectation p_i/p_j m_ij and a binomial distributed random variable with parameters n_i and 1/n_in_i/n_j  m_ij is O(1/√(n)). Hence, as long as the number of particles born in any of the two branching processes is o(√(n)), the two processes can be perfectly coupled w.h.p. In the remainder of this section we show that the coupling is w.h.p. perfect up to time t^* with t^* given by (<ref>). 



By construction, if the first |V̂^a(ℓ^*)| + |V̂^e(ℓ^*)| vertices that we “try to include” in  {V̂^a(ℓ) ∪V̂^e(ℓ);ℓ∈ℕ_0}  are all different and none of the first |V̂^p(ℓ^*)| that we “try to include” in  {V̂^p(ℓ);ℓ∈ℕ_0} are in  V̂^a(ℓ^*) ∪V̂^e(ℓ^*), then the process is not flagged. 
 


The law of large numbers yields 

    |V̂^a(ℓ^*)| ≤ 2 max_i,j ∈ [K]n_i/n_jm_ij  |V̂^e(ℓ^*)|,

Here we use that the expected number of edges in E' with any given vertex in V' as head is bounded from above by max_i,j ∈ [K]n_i/n_jm_ij. Equation (<ref>) implies that 

    |V̂^e(ℓ^*)| + |V̂^a(ℓ^*)|  = O(|V̂^e(ℓ^*)|),


Using Assumption <ref> we obtain

    |V̂^p(ℓ^*)| ≤ n^κ (|V̂^a(ℓ^*)|+ |V̂^e(ℓ^*)|),

Combining inequalities (<ref>) and (<ref>) yields

    |V̂^p(ℓ^*)| = O_p(n^κ|V̂^e(ℓ^*)|),


Since {𝒵(t);t ≥ 0} can be perfectly coupled with {𝒮^(n)_v(t);t ≥ 0} w.h.p. until o(√(n)) particles are born in {𝒵(t);t ≥ 0}, we obtain by using (<ref>),  that |𝒵(t^*)| = o_p(n^(1-κ)/3). Therefore, we also know that |𝒵(t^*)| = o(n^(1/2) and |V̂^e(ℓ^*)| = o_p(n^(1-κ)/3).
Combined with (<ref>) this gives

    |V̂^e(ℓ^*)| + |V̂^a(ℓ^*)| = o_p(n^(1-κ)/3)

and 

    |V̂^p(ℓ^*)| = o(n^(1+2κ)/3) = o_p(n).

Using birthday-problem-like arguments <cit.> the probability that the first |V̂^e(ℓ^*)| + |V̂^a(ℓ^*)| activated vertices in {Ĝ(ℓ); ℓ∈ℕ_0} are not all different is bounded from above by 

    (min_j ∈ [K]n_j)^-1(|V̂^e(ℓ^*)| + |V̂^a(ℓ^*)|)^2= o_p(n^-(1+2κ)/3),

while the probability that among the first |V̂^p(ℓ^*)| vertices that we “try to include” in  {V̂^p(ℓ);ℓ∈ℕ_0} there are vertices in  V̂^a(ℓ^*) ∪V̂^e(ℓ^*) is bounded from above by 

    |V̂^p(ℓ^*)| |V̂^e(ℓ^*)| + |V̂^a(ℓ^*)|/min_j ∈ [K]n_j   = o_p(n^(1+2κ)/3) × o_p(n^(1-κ)/3) × O(n^-1) 
       =o_p(n^(κ-1)/3)
       =o_p(1).

We conclude that the probability that the exploration process  {Ĝ(ℓ); ℓ∈ℕ_0} is flagged up to and including step ℓ^* goes to 0 as n →∞, as we desired. We summarise the main result of this section in the following lemma.


There exists a probability space on which we can define the branching process {𝒵(t);t ≥ 0} and the susceptibility process {𝒮^(n)_v(t);t ≥ 0} for all n, such that

    ℙ(|𝒮^(n)_v(t;a,j)| = Ẑ_j^i(t; a)  ) → 1,

where t^*= t^*(n) is given by (<ref>).





§ PROOFS




 §.§ Proof of Theorem <ref>

We are interested in the expected fraction of vertices infected by vertices of type i among the ultimately infected vertices of type j during a major outbreak. By exchangeability this expected fraction is given by 

    ℙ(|,v ∈ V_j).









Consider 𝒮^(n)_v(t^*), where we assign to each vertex v' ∈𝒮^(n)_v(t^*) the value σ'_v' = d(v',v). In Section <ref> we constructed part of the graph G' by looking backward in time. To analyse (<ref>), we construct another part of the graph G', by looking forward in time and starting at V_init. By Assumption <ref>, 𝒮^(n)_v(t^*) does not overlap with V_init w.h.p. We condition on this event.

We construct the relevant part of G' through a series of subgraphs {G̃(ℓ); ℓ∈ℕ_0} as follows.



  1) Construct 𝒮^(n)_v(t^*) and the edges that connect these vertices in G' as in Section <ref>. Let G̃(0) be this graph together with (the isolated) vertices in V_init. 


  2) The vertices in V_init are active in G̃(0). 


  3) Assume that we know G̃(ℓ). Let σ̃_v(ℓ) be the distance from V_init to v in G̃(ℓ). If there is no path from V_init to v in G̃(ℓ) then we set σ̃_v(ℓ) = ∞. If there are no active vertices in G̃(ℓ) with distance from V_init less than σ̃_v(ℓ)-t^*, then the shortest path from V_init to v is the same in G' as in G̃(ℓ), and we set G̃(k) = G̃(ℓ), for all k ≥ℓ.


  4) If there are active vertices in G̃(ℓ) with distance from V_init less than σ̃_v(ℓ)-t^*, then we construct G̃(ℓ+1) as follows. Pick the active vertex with the lowest distance from V_init in G̃(ℓ) (in case of a tie, which occurs if there are still several active vertices in V_init, make a uniform choice among those vertices). Say that this vertex is u ∈ V_i. Then assign to u the point processes (ξ_u^1, ⋯ξ_u^K), having the correct distribution (see Section <ref>) and label the points of ξ_u^j (j∈ [K]) with independent uniform vertices from V_j without replacement.

First we check whether some of the labels chosen correspond to vertices in 𝒮^(n)_v(t^*). If there are such labels, we add “preliminary” edges with tail u and heads equal to the respective vertices in 𝒮^(n)_v(t^*) (whether the preliminary edges become “actual” edges in the graph is determined in step 5 of the construction). The lengths of such edges correspond to the points in the point processes. 
Say that u' ∈𝒮^(n)_v(t^*) ∩ V_j is one of the labels chosen and that u is chosen at “age” t'. 

If u' ∉𝒮^(n)_v(t^*;t',j) for at least one of the vertices that u connects to, then the distance from u to v becomes less than t^* in the preliminary graph. We set G̃(ℓ+1) = G̃(ℓ). Note that  we have identified all vertices from which there is a path of length at most t^* to v in the exploration of 𝒮^(n)_v(t^*). However, since we know that u ∉𝒮^(n)_v(t^*), we have to condition on the event that no edges with tail u and heads in the set 𝒮^(n)_v(t^*-t') is shorter than t'.


  5) Finally, in deciding whether the preliminary edges become edges in the graph G̃(ℓ), we consider the following. If the distance from u to v does not become less than t^* in the preliminary graph, then u becomes passive in G̃(ℓ+1), and the edge set of G̃(ℓ+1) consists of all edges in G̃(ℓ) plus the edges with tail u and heads corresponding to the labels of the points in (ξ_u^1, ⋯ξ_u^K). The edges that are added have the obvious lengths. The heads of the edges in G̃(ℓ+1) are part of the vertex set of G̃(ℓ+1). The active vertices in G̃(ℓ+1) are the active vertices in G̃(ℓ) apart from u plus the newly added vertices in G̃(ℓ+1).




In the construction above, vertices of type j are added in an interchangeable way, j∈[K]. If a preliminary edge of length a has a head in 𝒮^(n)_v(t^*)∩ V_j, then every vertex in 𝒮^(n)_v(t^*;a,j) has the same probability of being the head of this edge and becoming part of the constructed graph.

Next, for the coupling, we analyse 𝒮^(n)_v(t^*;a,j) even further. Denote the vertices in {u ∈ V_i; (u,v) ∈ E'}, i.e. the vertices in V_i, which are in the first generation of the susceptibility set of v, by u_i1,u_i2, ⋯, u_i, X_i(v). Here 
    X_i(v) = |{u ∈ V_i; (u,v) ∈ E'}|
 is the number of vertices of type i in this first generation of the susceptibility set. Note that X_i(v) is distributed as X_ij where X_ij is
binomially distributed with parameters n_i and m_ij/p_j. So, X_ij converges in distribution to a
Poisson distributed random variable with expectation (p_i/p_j)m_ij. For r ∈ [X_i(v)], let τ_ir(v)= η(u_ir,v) be the length of edge (u_ir,v) ∈ E'. Here τ_ir(v) is distributed as τ_ijr with distribution function ℙ(τ_ijr≤ a)=𝔼[ξ_ij(a)]/m_ij. For convenience, we let X_i=X_i(v) and τ_ir=τ_ir(v).

We consider the susceptibility processes of the vertices in {u_ir; i ∈ [K], r ∈ [X_i]} up to distance t^* from v separately, i.e. we consider the susceptibility set of u_ir up to time t^*- τ_ir. Those susceptibility sets are (following the arguments in Section <ref>) w.h.p.  not overlapping and independent. Also note that

    𝒮^(n)_v(t^*;a,j)= v ∪⋃_i =1^K ⋃_r=1^X_i𝒮^(n)_u_ir(t^*-τ_ir;a,j).

By Lemma <ref>, we know that 

    |𝒮^(n)_u_ir(t^*-τ_ir,a,j)| = Ẑ^i_j(t^*-τ_ir; a)

Furthermore, by equations (<ref>) and  (<ref>), we obtain that 

    e^-α tẐ^i_j(t;  a) → c(a,j) W^i

on {Ẑ^i(t) →∞}. Since the sets 𝒮^(n)_u_ir(t^*-τ_ir) are w.h.p. not overlapping (and do not contain v), it follows from (<ref>) that

    |𝒮^(n)_v(t^*;a,j)∖ v|= ∑_i =1^K ∑_r=1^X_i |𝒮^(n)_u_ir(t^*-τ_ir;a,j)|

Then, by (<ref>) and (<ref>), for all i ∈ [K] and r ∈ [m_i],

    e^-α t^* |𝒮^(n)_u_ir(t^*-τ_ir;a,j)| → e^-ατ_ir  c(a,j) W^i(r)

in probability on  {Ẑ^i_j(t) →∞} as n →∞. Here, the random variables W^i(1), W^i(2), … are independent copies of W^i, i∈[K]. Let Ŵ^i be distributed as W^i(Ẑ^i(t) →∞), and let Ŵ^i(r) be defined analogously. Note that {Ẑ^i(t) →∞} implies that the branching process survives and new particles are born w.h.p. in the interval [t^*/2,t^*]. If {Ẑ^i(t) ↛∞}, then there is a last birth in the process. Because t^*/2 →∞ as n →∞, there is then no particle born in in the interval [t^*/2,t^*] w.h.p. It follows from the coupling arguments above that, as n →∞,

    ℙ(𝒮^(n)_v(t^*) ≠𝒮^(n)_v(t^*/2)) →ℙ(Ẑ^i(t) →∞).


We are interested in the fraction of vertices (possibly specified by type and age) at time t^* in 𝒮^(n)_v(t^*) that are connected to v through a path of vertices in 𝒮^(n)_v(t^*) that include u_ir, i ∈ [K] and r ∈ [X_ij]. That is, we want to analyse 

    |𝒮^(n)_u_ir(t^*-τ_1r;a,j)|/∑_j =1^K ∑_r=1^X_j |𝒮^(n)_u_jr(t^*-τ_jr;a,j)|,

on the set ∑_j =1^K ∑_r=1^X_j |𝒮^(n)_u_jr(t^*-τ_jr;a,j)| ≠ 0. Lemma <ref> allows us to couple the epidemic process with a branching processes {𝒵(t);t ≥ 0} such that w.h.p.

    |𝒮^(n)_u_ir(t^*-τ_1r;a,j)|/∑_j =1^K ∑_r=1^X_j |𝒮^(n)_u_jr(t^*-τ_jr;a,j)| = 
    Z_j^i,r(t^*-τ_1r;a)/∑_k =1^K ∑_r=1^X_kj Z_j^k,r(t^*-τ_kr;a),

where Z_j^i,r(t;a) are independent copies of Z_j^i(t;a), r ∈ℕ. Multiplying numerator and denominator of (<ref>) by e^-α t^* and using (<ref>) we obtain that (<ref>) is equal to

    e^-ατ_ir (e^-α (t^*-τ_ir)Z_j^i,r(t^*-τ_1r;a))/∑_k =1^K ∑_r=1^X_kj e^-ατ_kr (e^-α (t^*-τ_kr) Z_j^k,r(t^*-τ_kr;a),

which converges a.s. to

    e^-ατ_ir c(a,j) Ŵ^i(r)/∑_k =1^K ∑_r=1^X_kj e^-ατ_kr c(a,j) Ŵ^k(r) = e^-ατ_irŴ^i(r)/∑_k =1^K ∑_r=1^X_kj e^-ατ_krŴ^k(r)

as n →∞ (i.e. as t^* →∞). Note that the right hand side of (<ref>) does not depend on a or j. Since vertices of type j, j ∈ [K], are added in such a way that every vertex in 𝒮^(n)_v(t^*;a,j) has the same probability of being the head of this edge, every vertex in 𝒮^(n)_v(t^*;a,j) has the same probability of becoming part of the constructed graph. So, 
ℙ(|) is equal to

    1/ℙ(𝒮^(n)_v(t^*) ≠𝒮^(n)_v(t^*/2))𝔼[∑_r=1^X_i |𝒮^(n)_u_ir(t^*-τ_ir;a,j)|/∑_k =1^K ∑_r=1^X_k |𝒮^(n)_u_kr(t^*-τ_kr;a,j)|(𝒮^(n)_v(t^*) ≠𝒮^(n)_v(t^*/2))
    ]

and converges to ρ_ij, with ρ_ij given by (<ref>). This completes the proof for Theorem <ref>.












 §.§ Towards to proof of Theorem <ref>: Bounds for the ρ_ij

Theorem <ref> provides us with an expression for the asymptotic fractions ρ_ij of infected individuals of type j that were infected by individuals of type i. However, often there is no explicit description of the distribution of W^k(r). In Theorem <ref> we consider bounds for the ρ_ij for a special class of models, as specified in Section <ref>. In order to obtain those bounds, in this subsection and Section <ref>, we discuss, for the general setting, how to obtain the bounds for ρ_ij using the epidemic random graph G'. In most of the subsequent analysis we analyse the graph G' without taking the lengths of edges in into account. 







Note that Assumption <ref> guarantees that, with probability 1, all paths in E' have different lengths. In order to obtain the maximum and minimum of the probability ρ_ij for fixed mean offspring matrix M, we first investigate the susceptibility set of v_* (v_* ∈ V_j) restricted to the graph G'_ij=(V,E'_ij). We denote this susceptibility set by 𝒮_v_*,ij. Here E'_ij⊂ E' is the subset of E' that consists of the edges that either have tail vertex in V_i or head vertex not in V_j, i.e. E ∖ E'_ij is the set of edges with tails in V ∖ V_i and heads in V_j.
In Figure <ref>, the set 𝒮_v_*,11 consist of all  vertices of type 1 (the circles).

Let 𝒮^j_v_*,ij= V_j ∩𝒮_v_*,ij. If 𝒮_v_*,ij∩ V_init = ∅, i.e. if there is no path from V_init to v_* in G'_ij, then, by the definition of the epidemic process, every  vertex of type j in 𝒮_v_*,ij has the same probability to be the first one to be infected in the epidemic, i.e. u_*=_ u ∈𝒮^j_v_*,ijd(V_init,u) is uniform in 𝒮^j_v_*,ij.

Condition on 𝒮_v_*∩ V_init≠∅. If u_* =v_* and 𝒮_v_*,ij∩ V_init = ∅, then v_* is infected by a vertex that is not of type i. On the other hand, if u_* ≠ v_* or 𝒮_v_*,ij∩ V_init≠∅, then v_* might be infected by a type i vertex. Hence

    1-ρ_ij   = ℙ(|𝒮_v_*∩ V_init≠∅)
       ≥ℙ(u_* = v_*|𝒮_v_*∩ V_init≠∅).


Now assume that ℙ(η_i',j'> n^-1) =1, for i' ∈ [K]∖ i and j'=j, with j the type of v_*, and ℙ(η_i',j< n^-2) =1 otherwise. This implies that the lengths of any path with only edges in E'_ij is less than any edge in E'∖ E'_ij. The assumptions guarantee that, for u_* ≠ v_* or 𝒮_v_*,ij∩ V_init≠∅, conditioned on 𝒮_v_*∩ V_init≠∅, v_* is infected by a vertex of type i (the tail of the edge with head v_* in the shortest path in E_ij' from u_* or V_init to v_*). Therefore, for this model,

    ℙ(|𝒮_v_*∩ V_init≠∅) = ℙ(u_* = v_*|𝒮_v_*∩ V_init≠∅).

Hence, for a given distribution of E', models with ℙ(η_i',j'> n^-1) =1, for i' ∈ [K]∖ i and j'=j  and ℙ(η_i',j< n^-2) =1 otherwise, are among the models for which the fraction of the ultimately infected vertices of type j, infected by a type i vertex is maximised. 

Using similar arguments we obtain that ρ_ij, the fraction of ultimately infected  vertices of type j that are infected by vertices of type i, is minimal if the edge lengths of vertices with tail in i and head in j are much longer than the other edges. This will be used in Section <ref>.  





 §.§ Proof of Theorem <ref>

We consider the so-called symptom-response SEIR epidemic model that is introduced in detail in <cit.>, see also Remark <ref>.
That is, consider the model introduced in Section <ref> with K=2 and consider general distributions (ξ_i1, ξ_i2). 
For this model, using the arguments from Section <ref>, we can compute the maximal probabilities ρ_11^+, ρ_22^+ (and minimal probabilities ρ_21^- and ρ_11^-) explicitly. 
In this subsection we consider the model in which edges from V_1 to V_1 are infinitesimally short and all other edges in E' are relatively long. 
From Section <ref> we know that this is the model for which the fraction of vertices of type 1 infected by type 1 vertices is maximised. 
For reasons of convenience we assume that v_* is of type 1 with 𝒮_v_*∩ V_init≠∅. Note that Theorem <ref> considers the special case that (ξ_i1, ξ_i2) is obtained from the independent labelling of a one-dimensional point process ξ_i. In this special case we have ρ_11=ρ_1=ρ_12 and ρ_21=ρ_2=ρ_22. We treat Theorem <ref> in Remark <ref> at the end of this section.
First, we compute the upper bound ρ_11^+ (and lower bound ρ_21^-) for the general setting. 
We note that it is harder (if not impossible) to obtain an explicit expression for ρ_12^+ or ρ_21^+. As will become clear in the computation below, the difficulty with ρ_12^+ or ρ_21^+ is that one would need to consider paths in G_12' and G_21' that contain vertices of both type 1 and type 2. In contrast, for computing ρ_11^+,  paths in G_11' that end in v^* contain only vertices of  type 1.





If we ignore the lengths of edges in G' in the general model introduced in Section <ref>, then
the approximating (backward) branching process describing the generation-based growth of 𝒮_v_*(t) is defined through the following offspring distributions. The number of children of type i of a particle of type j is Poisson distributed with expectation m_ji^(b) = p_i/p_jm_ij, i,j∈[2]. For different i and j the distributions are independent of each other.

Let Y= |𝒮_v_*,11|. It is easily seen that Y is approximated by the size of a branching process with Poisson offspring distribution that has expectation m_11^(b). Then Y is Borel distributed with parameter m_11^(b), i.e. for ℓ∈{1,2,⋯}, 

    ℙ(Y=ℓ) = (m_11^(b)ℓ)^ℓ-1 e^-m_11^(b)ℓ/ℓ!

(see <cit.>). If m_11^(b)>1, then ℙ(Y=∞) > 0. Standard results on Borel distributions <cit.> give that, for m_11^(b)≤ 1,

    𝔼[1/Y] = 1-m_11^(b)/2.

Define

    ρ_11 = ℙ(|𝒮_v_*∩ V_init≠∅),

for v_* ∈ V_1. From the arguments in Section <ref> we know that

    ρ_21^- = 𝔼[1/Y|𝒮_v_*∩ V_init≠∅.]

Consistency then yields

    ρ_11^+=1-ρ_21^-.



We use the following extinction probabilities in the backward branching process:


  q_1: the probability that the backward branching process starting with a single type 1 particle goes extinct,

  q_2: the probability that the backward branching process starting with a single type 2 particle goes extinct,

  q̃_1: the probability that the backward branching process restricted to type 1 particles dies out, i.e. q̃_1 = ℙ(Y<∞).


From theory on multi-type supercritical branching processes <cit.> we know that (q_1, q_2) is the unique solution in (0,1)^2 of

    x     =    ∑_k=0^∞(m_11^(b))^k/k!e^-m_11^(b) x^k
    ∑_ℓ=0^∞(m_12^(b))^ℓ/ℓ!e^-m_12^(b) y^ℓ 
    = e^-[m_11^(b)(1-x)+ m_12^(b)(1-y)]
    
    y     =    ∑_k=0^∞(m_21^(b))^k/k!e^-m_21^(b) x^k
    ∑_ℓ=0^∞(m_22^(b))^ℓ/ℓ!e^-m_22^(b) y^ℓ 
    = e^-[m_21^(b)(1-x)+ m_22^(b)(1-y)].

Furthermore, q̃_1 is the smallest positive solution of (<ref>) which might be 1 or strictly smaller than 1, depending on whether or not p_1m_11 = m_11^(b)≤ 1.




 

We distinguish between m_11^(b)≤ 1 and m_11^(b)>1. This distinction is not necessary, but we think the argument becomes clearer, by treating the case q̃_1= 1 separately.




  §.§.§ The case m_11^(b)≤ 1

Assume that m_11^(b)≤ 1. Let 𝒜 be the event that the backward branching process with ancestor of type 1, involving both type 1 and type 2 individuals, survives. We explore the backward branching process of a particle of type 1 on 𝒜 as follows.




  * Explore the backward process of particles of type 1. If we ignore the conditioning on 𝒜, the process can be described by a subcritical branching process with Poisson offspring distribution with expectation m_11^(b). The random variable  Y is the total size of this branching process, including the initial individual. We know that Y is Borel(m_11^(b)) distributed.  



  * Condition on event 𝒜: the probability that a particle of type 1 has infinitely many descendants is ℙ(𝒜)= 1-q_1, where q_1 is defined through (<ref>) and (<ref>). We also use the probability that a particle has infinite offspring, conditioned on having no children of type 1. This probability is 1-q, where

    q= ∑_ℓ =0^∞(m_12^(b))^ℓ/ℓ ! e^-m_12^(b) (q_2)^ℓ = e^-m_12^(b)(1-q_2)=q_1 e^m_11^(b) (1-q_1).

Using Bayes' rule we obtain

    ℙ(Y=ℓ|𝒜)     = ℙ(𝒜|Y=ℓ) ℙ(Y= ℓ)/1-q_1
        = 1-q^ℓ/1-q_1(m_11^(bℓ)^ℓ-1 e^-m_11^(b)ℓ/ℓ!
       = 1/1-q_1(m_11^(b)ℓ)^ℓ-1 e^-m_11^(b)ℓ/ℓ! - q/1-q_1(m_11^(b)ℓ q)^ℓ-1 e^-m_11^(b)ℓ/ℓ!
       =1/1-q_1(  (m_11^(b)ℓ)^ℓ-1 e^-m_11^(b)ℓ/ℓ! -  q_1(m_11^(b) q_1ℓ)^ℓ-1 e^-m_11^(b) q_1 ℓ/ℓ!).

where we used (<ref>) in the last equality.


  * Using (<ref>), (<ref>)  and (<ref>) we find that 

    ρ_21^-   = 𝔼[Y^-1|𝒜]
       = 1/1-q_1(1- m_11^(b)/2 -q_1(1-m_11^(b) q_1/2) )
       = 1/1-q_1(1- q_1 - (1-(q_1)^2) m_11^(b)/2)
       = 1 - (1+q_1) m_11^(b)/2.



  * Using (<ref>) we find the desired expression ρ_11^+ for m_11^(b):

    ρ_11^+=1-(1+q_1) m_11^(b)/2.







  §.§.§ The case m_11^(b)  > 1

Assume that m_11^(b)  > 1. Then it is possible that the backward process restricted to particles of type 1 is already large, i.e. the approximating branching process with Poisson(m_11^(b)) offspring distribution already survives (call this event 𝒜_1). In that case the probability that vertex v is infected by a vertex of type 1 approaches 1 as the population size tends to infinity. 

The other possibility is that the backward process restricted to particles of type 1 stays small. Call this event 𝒜_1^C, the complement of 𝒜_1. The probability of this event is ℙ(𝒜_1^C)=q̃_1, where q̃_1 is the unique solution in (0,1) of equation (<ref>) (note that m_11^(b)>1 ensures that q̃_1 exists). From the theory of branching processes we know that, conditioned on 𝒜_1^C, the approximating branching process is still a branching process with Poisson distributed offspring distribution, but now with offspring expectation m_11^(b)q̃_1. We still condition on the event that the approximating branching process (including both types of particles) survives.

Conditioned on 𝒜_1^C, the total size of the backward process restricted to particles of type 1 is Borel distributed with parameter m_11^(b)q̃_1. We use Bayes' rule, 

    ℙ(Y=ℓ|𝒜,𝒜_1^C) = ℙ(𝒜|Y=ℓ,𝒜_1^C) ℙ(Y= ℓ|𝒜_1^C)/ℙ(𝒜|𝒜_1^C).

Note that 

    ℙ(𝒜|𝒜_1^C) = ℙ(𝒜,𝒜_1^C)/ℙ(𝒜_1^C)= ℙ(𝒜_1^C)-ℙ(𝒜^C,𝒜_1^C)/ℙ(𝒜_1^C)= ℙ(𝒜_1^C)-ℙ(𝒜^C)/ℙ(𝒜_1^C)= q̃_1-q_1/q̃_1.

Then, with the same arguments as those leading to (<ref>) but with m_11^(b) replaced by m_11^(b)q̃_1, yield

    ℙ(𝒜|Y=ℓ,𝒜_S^C) =1- q^ℓ = 1-(q_1)^ℓ e^m_11^(b) (1-q_1)ℓ.

and 

    ℙ(Y= ℓ|𝒜_1^C)= 
    (m_11^(b)q̃_1 ℓ)^ℓ-1 e^-m_11^(b)q̃_1 ℓ/ℓ!.

Combining these identities with (<ref>) yields

    ℙ(Y=ℓ|𝒜,𝒜_1^C)
       = q̃_1/q̃_1-q_1(1-(q_1)^ℓ e^m_11^(b) (1-q_1)ℓ)((m_11^(b)q̃_1 ℓ)^ℓ-1 e^-m_11^(b)q̃_1 ℓ/ℓ!)
       =  q̃_1/q̃_1-q_1((m_11^(b)q̃_1 ℓ)^ℓ-1 e^-m_11^(b)q̃_1 ℓ/ℓ!) - q_1 (q̃_1)^ℓ/q̃_1-q_1((m_11^(b) q_1 ℓ)^ℓ-1 e^-m_11^(b) (q̃_1+q_1-1) ℓ/ℓ!)
        =  q̃_1/q̃_1-q_1((m_11^(b)q̃_1 ℓ)^ℓ-1 e^-(m_11^(b)q̃_1 ℓ/ℓ!) - q̃_1/q̃_1-q_1((m_11^(b) q_1 ℓ)^ℓ-1 e^-m_11^(b) q_1 ℓ/ℓ!).

By standard results on Galton Watson branching processes we know that m_11^(b)q̃_1 ≤ 1. Moreover, since q_1 ≤q̃_1, also m_11^(b) q_1 ≤ 1 holds. Then, using (<ref>), we obtain

    𝔼[Y^-1|𝒜,𝒜_1^C]    = q̃_1/q̃_1-q_1(1-m_11^(b)q̃_1/2) - 
     q_1/q̃_1-q_1(1-m_11^(b) q_1/2)
        = 1 -  m_11^(b) (q̃_1 + q_1)/2.

This leads to 

    ρ_21^-    = 𝔼[Y^-1|𝒜]
       = 𝔼[Y^-1|𝒜,𝒜_1^C]ℙ[𝒜_1^C|𝒜] + 𝔼[Y^-1|𝒜,𝒜_1] ℙ[𝒜_1|𝒜]
       = 𝔼[Y^-1|𝒜,𝒜_1^C]ℙ[𝒜_1^C,𝒜]/ℙ[𝒜] + 0 
       = 𝔼[Y^-1|𝒜,𝒜_1^C]ℙ[𝒜_1^C]-ℙ[𝒜_1^C,𝒜^C]/ℙ[𝒜]
       = ( 1 -  m_11^(b) (q̃_1 + q_1)/2) q̃_1-q_1/1-q_1.

Note that this expression (<ref>) is consistent with the result (<ref>) for m_11^(b)≤ 1, where q̃_1=1. Finally, using (<ref>), we find the expression for ρ_11^+:

    ρ_11^+=1-( 1 - m_11^(b) (q̃_1 + q_1)/2) q̃_1-q_1/1-q_1.






In <cit.> we interpret vertices of type 1 as individuals that show symptoms when infectious, while vertices of type 2 are asymptomatic throughout their infectious period. We assume that whether infected individuals become symptomatic or not, does not depend on who infected them or when they were infected. So, we may assign i.i.d. types to the vertices before the epidemic. Note that the type of a vertex that does not become infected during the epidemic has no epidemiological relevance. Type 1 and type 2 vertices are equally susceptible. This implies that {(ξ_v^1(t),ξ_v^2(t)); t ≥ 0} can be obtained by considering a one-dimensional point process {ξ_v(t); t ≥ 0}, for all v∈ V.  One can assign the types (type 1 with probability p_1 or type 2 with probability p_2=1-p_1) independently to the points of this process. Then ρ_11=ρ_1=ρ_12 and ρ_22=ρ_2=ρ_21. Furthermore m_ji^(b)=p_i m̃_j, with m̃_j=𝔼[ξ_j(∞)] the expected number of secondary cases generated by a newly infected type j individual in an otherwise susceptible population, i,j=1,2. Note that (<ref>) reduces to the fraction ρ_1^+ of the ultimately infected vertices, that are infected by symptomatic vertices: 

    ρ_1^+ = 1-( 1 - p_1m̃_1(q̃_1 + q)/2) q̃_1-q/1-q,

where q̃_1 and q are solutions of (<ref>) and (<ref>). This proves Theorem <ref>.





§ DISCUSSION


In this manuscript we couple a (fairly) general multi-type stochastic epidemic process with a weighted random graph. We use this random graph to obtain a characterisation of the (large population limit) fraction of individuals in the population that had an infector of (say) type 1 given a large outbreak. 
The results of this paper are applied (and in more detail in <cit.>) to a model where the types of individuals represent whether the infected individual will show symptoms at some moment if infected, or he or she stays asymptomatic  if  infected. 

From a public health perspective, a relevant variant of the model would be to consider a population with only one type of individuals, where infectious individuals may start off asymptomatic after which they become symptomatic (see <cit.>, where the asymptomatic phase is also referred to as “presymptomatic”). The related question for this model would be: “what fraction of the infected population had a symptomatic infector, given a large outbreak?”. In order to analyse this model we need to assign types to the edges in the epidemic graph instead of to the vertices. We can use the techniques of this paper to characterise the answer to the above question. We can still define a susceptibility process from a vertex v. The susceptibility process can be approximated by a multi-type branching process. The type of a particle in the branching process should then corresponds to the type of the edge through which the particle is added to the susceptibility process. That is, the type of the particle (say u) in the branching process depends on whether the first edge in the shortest path from u to v in the epidemic graph  represents a contact that was made while u was symptomatic or asymptomatic. As in Theorem <ref> the answer will be implicit through its dependence on the distribution of the martingale limits W.



§ ACKNOWLEDGEMENTS

The authors are supported by Vetenskapsrådet (Swedish Research Council), grants 2015-05015
(TB and KYL) and  2016-04566 (PT).

 

abbrv


