Modeling the early transmission of COVID-19 in New York and San Francisco using a pairwise network model.

Classical epidemiological models assume mass action. However, this assumption is violated when interactions are not random. With the recent COVID-19 pandemic, and resulting shelter in place social distancing directives, mass action models must be modified to account for limited social interactions. In this paper we apply a pairwise network model with moment closure to study the early transmission of COVID-19 in New York and San Francisco and to investigate the factors determining the severity and duration of outbreak in these two cities. In particular, we consider the role of population density, transmission rates and social distancing on the disease dynamics and outcomes. Sensitivity analysis shows that there is a strongly negative correlation between the clustering coefficient in the pairwise model and the basic reproduction number and the effective reproduction number. The shelter in place policy makes the clustering coefficient increase thereby reducing the basic reproduction number and the effective reproduction number. By switching population densities in New York and San Francisco we demonstrate how the outbreak would progress if New York had the same density as San Francisco and vice-versa. The results underscore the crucial role that population density has in the epidemic outcomes. We also show that under the assumption of no further changes in policy or transmission dynamics not lifting the shelter in place policy would have little effect on final outbreak size in New York, but would reduce the final size in San Francisco by 97%.

Introduction

The number of infected cases and deaths from the coronavirus disease 2019 (COVID-19) continues to increase (https://www.who.int/.; Li et al., 2020a). Since the first case of COVID-19 was confirmed in the United States on January 22th, 2020, the epidemic has spread in North America (https://www.who.int/.). There have been more than 16,000 daily new cases almost every day in the United States since March 23rd, 2020, and the number of total confirmed cases has continued to grow. In this paper, we investigate the early COVID-19 infections in the United States, explore the factors that are important in determining disease outcomes. Our approach is to consider two major US cities with very different infection histories: New York and San Francisco. We compare and contrast the dynamics in these two cities, with a focus on determining the demographic and behavioral factors driving the disease outbreak.

Coronavirus is transmitted mainly through close contact and respiratory droplets, with possible airborne transmission if aerosols are present (Chan et al., 2020; Huang et al., 2020). Available evidence shows that the transmission can arise from symptomatic, pre-symptomatic (exposed) and asymptomatic people infected with COVID-19 (Kimball et al., 2020; Ong et al., 2020; Pan et al., 2020; Yu et al., 2020). Symptoms include cough, fever and difficulty breathing, as with SARS (Severe acute respiratory syndrome coronavirus) and MERS (Middle East respiratory syndrome coronavirus) (Wu et al., 2020). The incubation period for COVID-19, is on average 5–6 days, although it can take up to 14 days. During this “pre-symptomatic” period, some infected people can be contagious. Therefore, transmission from pre-symptomatic cases can occur before symptom onset. It is also possible for infected people to remain asymptomatic while still transmitting the disease.

With the advent of COVID-19 there has been an explosion in the number of mathematical models to study its dynamics. Many of the early models involve analysis of the COVID-19 outbreak in China. A deterministic compartmental model, based on the clinical progression of the disease, epidemiological status of the individuals, and intervention measures to estimate the transmission risk of COVID-19 at the early stage of epidemic was proposed and analyzed by Tang et al. (Tang et al., 2020a, 2020b). Compartmental models have also been augmented with artificial intelligence (AI) approaches to predict the epidemics trend of COVID-19 in China (Yang et al., 2020).

As with other diseases, quarantining plays a major role in the COVID-19 outcomes. Hu et al. modelled quarantine strategies for COVID-19 at different population densities to show how this factor affects the disease outbreak in Guangdong, China (Hu et al., 2020). The effectiveness of quarantine and isolation in China was also discussed in some papers (Hou et al., 2020; Tang et al., 2020c).

One of the main policy strategies for addressing high COVID-19 infection levels is to use social distancing. The effects of such policies have been evaluated by mathematical models. Leung et al. used a susceptible-infectious-susceptible (SIS) model to show the potential effects of relaxing containment measures after the first wave of infection, in anticipation of a possible second wave (Leung et al., 2020). The impact of mixed social distancing and school closure on disease transmission was studied by Zhang et al. and it was found that social distancing alone is sufficient to control COVID-19 in China (Zhang et al., 2020). In other jurisdictions, Anderson et al. estimated the impact of physical distancing on the disease contact rate in British Columbia, Canada (Anderson et al., 2020).

One of the early models that most impacted US policy has come from Imperial College, London (Ferguson et al., 2020). There an individual-based model was used to explore scenarios for COVID-19 in Great Britain and the US, and it was shown that the health care system would be easily overwhelmed without sufficient suppression of the disease via social distancing, home isolation and household quarantine. To fight the COVID-19 pandemic, the United States has adopted a series of the above control measures. These including restricting travel and quarantining close contacts of those with COVID-19, as well as social distancing via advising people to stay at home and by keeping more than 2 m distance from others (https://www.cdc.gov/coron, 2019; https://news.sina.com.cn/, 2020; http://yl.szhk.com/2020/0, 2020).

Among these mathematical models for studying the transmission dynamics of COVID-19, the division of compartments varies. The simplest is to classify individuals as susceptible, infectious or recovered (Roda et al., 2020). It is more realistic to classify individuals as susceptible, exposed, infectious, and recovered. Considering that asymptomatic people infected with COVID-19 can transmit the disease, the infectious are further divided two types: the asymptomatic and the symptomatic (Tang et al., 2020a, 2020b). Given the policy of quarantining the close contacts with the infectious, people in different disease states are further divided into the unquarantined and the quarantined (Tang et al., 2020a, 2020b). To forecast the transmission of COVID-19 in the United States, it is more realistic to classify individuals as susceptible, exposed, asymptomatic infectious, symptomatic infectious, and recovered with considering of quarantine. Except for the compartments, a key factor in the transmission process is the number of contacts, since a susceptible person is infected by an infectious person at the same transmission rate per unit time in the same place. In a well-mixed compartment model, the number of contacts is constant (Roda et al., 2020), which is ideal and can not reflect the clustering effect of population. Although some compartment models considered that the number of contacts was time-varying (Tang et al., 2020a), it is assumed and man-made. The broken links caused by quarantine are not clear in these compartment models. In a social network, nodes represent individuals and links between nodes represent the contacts between two connected individuals. Network models consider the heterogeneity of contact structures, and use the degrees of nodes to reflect the number of contacts. Xue et al. used a degree-based network model to study the transmission of COVID-19 in Wuhan, Toronto and Italy (Xue et al., 2020). However, due to the policies of quarantine and reopening, the distribution of degree should be time-varying. The pairwise model regards the number of links as variables, and simultaneously studies the dynamic behavior of nodes and links between nodes, which effectively reflects human social behavior, and has natural advantages in studying quarantine and social distance. Quarantining the close contacts and keeping social distance lead to a decrease in the number of links, while resumption of work and school leads to an increase in the number of links. In fact, epidemiological surveys of close contacts of infected persons are based on their social networks. Compared with simple SIR or SEIR models, our pairwise network model is more likely to follow the realistic scenarios during the COVID-19 epidemic and to provide better quantitative results. By explicitly incorporating the effects of family cluster and contact tracing followed by household quarantine and isolation, our model provides a good fit to the trajectory of COVID-19 infections and rational prediction for the epidemic trend.

Here we focus our study on two different cities in the United States: New York and San Francisco. These two cities are chosen based on the following considerations. New York and San Francisco are two large United States cities, one initially hard-hit, the other initially not so hard hit. By looking at the differences between the two large cities we can try to understand factors that are important in determining disease outcomes. New York city includes 5 core counties: New York county (Manhattan), Kings county (Brooklyn), Bronx county (the Bronx), Richmond county (Staten Island), and Queens county (Queens). San Francisco city includes the greater area with 5 core counties: San Francisco, Alameda, Marin, Contra Costa, and San Mateo County.

In this work we build a pairwise network model, based on methods in Ref. (Luo et al., 2021), to investigate the spread of COVID-19 in New York and San Francisco. We use data of COVID-19 infections collected from Ref. (https://covidmapper.ca/ge) to identify the parameters of the proposed model and calculate some common descriptors of epidemic timing and impact, including the basic reproduction number and the effective reproduction number. The sensitivity analysis shows that there is a highly negative correlation between the clustering coefficient and the basic reproduction number. Comparing the transmission in New York and San Francisco, we find that population density has a major effect on disease transmission. When the population density of New York drops to the same as San Francisco, the epidemic duration will be shortened by more than 140 days, and the final size will be reduced by more than 50%. When the population density of San Francisco increases to the same as New York, the final size will be increased 15-fold. Not lifting the shelter in place policy has a little effect on final size in New York, but will reduce the final size in San Francisco by 97%.

Materials and methods

Data collection and description

We collect all the daily reported data on COVID-19 infections, including the laboratory confirmed cases, and death-caused by COVID-19 in two typical cities in the United States: New York city and San Francisco city (https://covidmapper.ca/ge). The data is from March 27th, 2020 to July 20th, 2020. We also collect population sizes and areas of these two cities (https://en.wikipedia.org/a; https://en.wikipedia.org/b; https://en.wikipedia.org/c; https://en.wikipedia.org/d; https://en.wikipedia.org/e; https://en.wikipedia.org/f). The population and area of New York are 8,336,697 and 785 square kilometers, respectively and those of San Francisco are 4,264,934 and 6405.5 square kilometers. This means that the population density of New York is about 16 times as dense as that of San Francisco.

In our social network, node represents individual and one link between two individuals describes that such two individuals spend more than 10 minutes at a distance of less than 2 m (https://www.cdc.gov/coron, 2019). The clustering coefficient denoted by φ, one of the network characteristic parameters, can characterize the isolation intensity of “stay at home” order (also called a shelter in place policy) (Luo et al., 2021). New York had a shelter in place policy from March 22nd to May 28th (https://rt.live/us/), and San Francisco had a shelter in place policy from March 19th to May 25th (https://rt.live/us/). The more intensive the household quarantine measure is, the larger the clustering coefficient is. According to Ref. (https://www.census.gov/da), the number of family members varies from 1 person to more than 7 persons in the United States. Since there is no data available on family structure in New York and San Francisco, we use the distribution of household numbers in the United States to describe the distribution in each of the two cities. The ideal case is that every family is completely isolated. In this case, every family is a complete graph.

Mathematically the clustering coefficient φ of a network is defined as the number of triangles over the number of triples in the network (Keeling, 1999). Therefore, in this ideal case, φ behaves aswhere F is the number of families each with k family members, C represents combination and A represents permutation. When the number of links with the outside for every family increases by three, ignoring the few addition of triangles, we obtain

We consider this case in the following.

The pairwise model for COVID-19 epidemic

To study the COVID-19 epidemic spread in the United States, based on Luo et al. (Luo et al., 2021), we propose a pairwise epidemic model with family clusters and quarantine of individuals having been contacted with the confirmed cases. As demonstrated in Fig. 1, individuals are grouped into 12 states: unquarantined susceptible (S), quarantined susceptible (S0), unquarantined incubation (E), quarantined incubation (E0), unquarantined asymptomatic (A), quarantined asymptomatic (A0), confirmed symptomatic (I), hospitalized (H), recovered asymptomatic (R), quarantined recovered asymptomatic (), recovered confirmed symptomatic (R) and died while recovering (R). The individuals in states of E, A, I and H are infectious, but H individuals do not transmit the disease as they have been isolated in hospital. The transitions between different states are listed as follows. Each S individual moves into the E state when it is exposed to the disease at the transmission rate β1 by an E individual, at rate β2 by an A individual and at rate β3 by an I individual. After an average incubation period of 1/q days, each E individual moves to I state with probability p or to A state with probability 1 − p. The individuals in state A enter into the R compartment at rate γ2, while the individuals in state I move into either the H compartment with proportion ϵ after 1/γ4 days or the R compartment with proportion 1 − ϵ after 1/γ1 days. When the I individuals are in the R compartment, they will reconnect k links on average with the individuals in states S, E, A, R or R. After an average period of 1/γ5 days, the H individuals move into R compartment with proportion 1 − d or R compartment with proportion d. When the H individuals move into R compartment, they will reconnect k links on average with the individuals in state S, E, A, R or R. Meanwhile, the individuals in states S, E, A or R will be quarantined with probability η and move to S0, E0, A0 or , once any one of their neighbors is confirmed to be infected. Accordingly, the individuals in state S0, E0, A0 or will go back to the S, E, A or R compartment after the quarantine period of average length 1/γ3 days. At the same time, they will reconnect k links on average with the individuals in states S, E, A, R or R. Since the individuals in state I are confirmed as symptomatic, we ignore the reconnection links with the individuals in this state.

Fig. 1

The flow diagram of pairwise model for COVID-19 epidemic in the United States.

Social distancing measures, such as keeping more than 2 m distance between individuals in public places, have contributed to the prevention and control of the COVID-19 epidemic. The average degree in our social networks describes the contact number of each individual with other individuals on average within 2 m distance, which is used to characterize the social distance among individuals in this paper. Keeping a social distance means breaking the links. Returning to work and school means reconnecting the links. For the change of social distance, there are two key elements: the intervention time of changing social distance among individuals and the intensity of intervention. We define an intervention function of social distance (IFSD):where δ(⋅) is the Kronecker delta function. Its numerical implementation is given as follows:where t > 0 is the intervention time of changing social distance among individuals, and α is the intensity of intervention. The time unit is day according to the daily data. Solving a system of differential equations generally requires discretization, therefore we consider that social distance changes during the interval [t, t + 1) with the intensity of intervention α ≠ 0 during this interval. After t + 1, the new social distance remains the same. It is worth noting that this function is defined in cases where social distance measures change quickly. Then the setting of intervention time makes sense. Take the resumption of work and school as an example. Once the time is set, almost all enterprises and schools open at the same time, and social networks have changed a lot at the same time. Even if some enterprises or schools are delayed for a few days, it will have little impact on the overall social network. The function θ(t; t) mirrors the change of social distancing measures. When α < 0 (the stricter social distancing measures), more links are broken, and thus the earlier the intervention time t the more effective COVID-19 is contained. In contrast, when α > 0 (relaxing social distancing measures such as work resumption), the later the intervention time t the more beneficial to curb the disease spread. Early resumption of work and school could easily lead to a second outbreak, whereas late resumption would do substantial damage to the economy. Therefore, the appropriate time and intensity of intervention are extremely important.

This pairwise model is applicable to a shelter in place policy. In New York, shelter in place started on March 22nd, and ended on May 28th. In San Francisco, shelter in place started on March 19th and ended on May 25th. When shelter ended, the resumption of work and school began. In this case, the contacts between individuals increase, and links in the social network were correspondingly reconnected, which leads to the larger average degree. In this case, the intensity of intervention α is positive.

Using the notation in Table 1 and parameters in Table 2, and considering the change of social distance measures, the pairwise model for COVID-19 epidemic can be described as follows:

Table 1

Model variables and their definitions.

variable	definition
[S]	The number of unquarantined susceptible individuals
[S0]	The number of quarantined susceptible individuals
[E]	The number of unquarantined exposed (incubation) individuals
[E0]	The number of quarantined exposed (incubation) individuals
[A]	The number of unquarantined asymptomatic individuals
[A0]	The number of quarantined asymptomatic individuals
[I]	The number of the confirmed symptomatic individuals
[H]	The number of the hospitalized individuals
[RI]	The number of recovered individuals with symptoms
[RA]	The number of recovered individuals without symptoms
[R0A]	The number of the quarantined recovered individuals without symptoms
[Rd]	The number of individuals who die while recovering
[SS]	Twice the number of links between nodes with S state and S state
[SE]	The number of links between nodes with S state and E state
[SA]	The number of links between nodes with S state and A state
[SI]	The number of links between nodes with S state and I state
[SRA]	The number of links between nodes with S state and RA state
[SRI]	The number of links between nodes with S state and RI state
[EE]	Twice the number of links between nodes with E state and E state
[EA]	The number of links between nodes with E state and A state
[EI]	The number of links between nodes with E state and I state
[ERA]	The number of links between nodes with E state and RA state
[ERI]	The number of links between nodes with E state and RI state
[AA]	Twice the number of links between nodes with A state and A state
[AI]	The number of links between nodes with A state and I state
[ARA]	The number of links between nodes with A state and RA state
[ARI]	The number of links between nodes with A state and RI state
[II]	Twice the number of links between nodes with I state and I state
[IRA]	The number of links between nodes with I state and RA state
[IRI]	The number of links between nodes with I state and RI state
[RARA]	Twice the number of links between nodes with RA state and RA state
[RARI]	The number of links between nodes with RA state and RI state
[RIRI]	Twice the number of links between nodes with RI state and RI state
[SSE]	The number of triples with the joint structure S–S-E
[SSA]	The number of triples with the joint structure S–S-A
[SSI]	The number of triples with the joint structure S–S–I
[SEE]	The number of triples with the joint structure S-E-E
[SAE]	The number of triples with the joint structure S-A-E
[SIE]	The number of triples with the joint structure S-I-E
[SRAE]	The number of triples with the joint structure S-RA-E
[ESE]	Twice the number of triples with the joint structure E-S-E
variable	definition
[ESA]	The number of triples with the joint structure E-S-A and [ASE] = [ESA]
[ESI]	The number of triples with the joint structure E-S-A and [ISE] = [ESI]
[ESRA]	The number of triples with the joint structure E-S-RA
[ESRI]	The number of triples with the joint structure E-S-RI
[EEE]	Twice the number of triples with joint structures E-E-E
[EEA]	The number of triples with the joint structure E-E-A
[EEI]	The number of triples with the joint structure E-E-I
[EERA]	The number of triples with joint structures E-E-RA
[EERI]	The number of triples with joint structures E-E-RI
[EAE]	Twice the number of triples with the joint structure E-A-E
[EARA]	The number of triples with the joint structure E-A-RA
[EARI]	The number of triples with the joint structure E-A-RI
[EIE]	Twice the number of triples with the joint structure E-I-E
[EIRA]	The number of triples with the joint structure E-I-RA
[EIRI]	The number of triples with the joint structure E-I-RI
[ERAE]	Twice the number of triples with joint structures E-RA-E
[ERARA]	Twice the number of triples with joint structures E-RA-RA
[ERARI]	Twice the number of triples with joint structures E-RA-RI
[ASA]	Twice the number of triples with the joint structure A-S-A
[ASI]	The number of triples with the joint structure A-S-I and [ASI] = [ISA]
[ASRA]	The number of triples with the joint structure A-S-RA
[ASRI]	The number of triples with the joint structure A-S-RI
[AAE]	The number of triples with the joint structure A-A-E
[AIE]	The number of triples with the joint structure A-I-E
[ARAE]	The number of triples with joint structures A-RA-E
[ISI]	Twice the number of triples with the joint structure I–S–I
[ISRA]	The number of triples with the joint structure I-S-RA
[ISRI]	The number of triples with the joint structure I-S-RI
[IIE]	The number of triples with the joint structure I–I-E
[IRAE]	The number of triples with the joint structure I-RA-E

Table 2

Definition of parameters and their values (unit time: day).

Parameter	Definition	New	Standard	San	Standard	Data
		York	deviation	Francisco	deviation	source
β1	transmission rate by E individuals	0.142 6	0.009 0	0.091 5	0.003 9	MCMC
β2	transmission rate by A individuals	0.228 4	0.014 3	0.229 9	0.005 7	MCMC
β3	transmission rate by I individuals	0.456 8	–	0.459 8	–	(Li et al., 2020b; Moghadas et al., 2020)
p	probability of showing symptoms	0.25	–	0.25	–	(https://www.cidrap.umn.ed, 2020; https://www.washingtonpos, 2020)
1/q	incubation period for E individuals	5	–	6	–	(Moss et al., 2020; Park et al., 2020)
η	quarantined probability of individuals who have been contacted with the confirmed cases	0.8	–	0.8	–	Moss et al. (2020)
ϵ	proportion of the confirmed symptomatic cases transferring to hospitalized cases	0.2	–	0.2	–	Moss et al. (2020)
1/γ1	period of transferring from I individuals directly to RI individuals	9.68	–	9.68	–	Moss et al. (2020)
γ2	The recovery rate of A individuals	1/12	–	1/12	–	Luo et al. (2021)
1/γ3	quarantine period for S, E, A and RA individuals	14	–	14	–	Moss et al. (2020)
1/γ4	period of transferring from I individuals directly to H individuals	4	–	4	–	Moss et al. (2020)
1/γ5	period of transferring from H individuals directly to RI individuals	5.68	–	5.68	–	Moss et al. (2020)
d	proportion of H individuals who die while recovering	0.075	–	0.075	–	(https://covidmapper.ca/ge)
N	size of total population	8,336,697	–	4,264,934	–	(https://en.wikipedia.org/a; https://en.wikipedia.org/b; https://en.wikipedia.org/c; https://en.wikipedia.org/d; https://en.wikipedia.org/e; https://en.wikipedia.org/f)
φ	clustering coefficient	0.297 5	–	0.297 5	–	(2)
k	the average degree of individuals released from quarantine	1.471 3	0.023 5	1.510 9	0.004 7	MCMC
α	the intensity of intervention	0.069 6	0.010 6	0.122 1	0.005 4	MCMC

Since the evolution of the number of individuals in each class with time in Eqs.(3), (4), (5), (6), (7), (8), (9), (10), (11), (12), (13), (14) depends on the numbers of different types of links, in order to close Eqs.(3), (4), (5), (6), (7), (8), (9), (10), (11), (12), (13), (14), we need to derive the evolution equations of the number of links of each type with time. Because individuals in states are quarantined, isolated or removed, there are no links among them. Let , the evolution equations of are as follows:

Here M = [S] + [E] + [A] + [R] + [R]. The terms represent the case as that individuals in state E turn into I, their neighbors are quarantined, resulting in a decrease in the number of links. The terms represent a lifting of the quarantine or that recovering results in an increase in the number of links. The change of social distance measures, reflected by the IFSD function, directly leads to the change of links in the social network. The IFSD function reflects both the intervention time t of changing social distance among individuals and the intensity of intervention α. The change of links is positively correlated with α, so we consider considered the simplest form θ(t; t)[XY] to represent the change of social distance measures. Since the individuals in state I are confirmed symptomatic, we ignore the change of links with the individuals in this state due to the change of social distance measures.

The numbers of triples could be closed by the following moment closure approximation formula (referred to in (House & Keeling, 2011; Keeling, Rand, & Morris, 1997)),where φ is the clustering coefficient, N is the number of nodes in the underlying network and n(t) is the average degree at time t with the initial condition n0 ≜ n(0), satisfying

As a pair approximation model, our model is somewhat different from that of Luo et al. (Luo et al., 2021). The main differences are listed as follows:

Our model considers the links reconnection between individuals released from quarantine, focusing on studying the effect of social distancing measures and work resumption on COVID-19 spread.

The quarantined probability of close contacts η is less than one, since close contacts with confirmed cases are quarantined at home rather than centrally.

Only about 20% are hospitalized (Moss et al., 2020), and most diagnosed patients stay at home and some of them may go out and infect others.

Our model considers more detailed recovered or removed compartments R, , R, and R.

The basic reproduction number and effective reproduction number

The basic reproduction number , which is defined as an average number of secondary cases that one case produces during the course of its infectious period in a totally susceptible population, determines the growing rate of an emerging infectious disease. Sometimes it can be used to predict whether an outbreak of one infectious disease will occur or not. Adopting the method proposed by the next generation matrix method (van den Driessche & Watmough, 2002), the expression takes the form ofwhereand

As almost all individuals are susceptible when the outbreak starts with the introduction of a single “average” infected individual at time t = 0, there is the constraint [SS](0) = n0[S](0) in networks.

Since the pairwise epidemic model (3–35) is a network model, we are mainly concerned with the variants of the infected links. The biological meaning of can be regarded as a sum of three quantities. R[ is the number of secondary S − E links that one S − E link will generate in an entirely S − S link group (1 − φ)(n0 − 1) during its lifespan . Note that and denote the probabilities that one S − E link survives to progress into S − A link and S − I link, respectively. Hence, R[ represents the number of secondary S − A links that one S − A link will produce in a totally S − S link group during its lifespan . R[ is the number of secondary S − I links that one S − I link will produce in a totally S − S link group during its lifespan .

The effective reproduction number is defined as an average number of secondary cases that one case produces during the course of its infectious period. Unlike the basic reproduction number, the effective reproduction number is given in real-time and depends upon the current state of the epidemic. Using the same method as with the basic reproduction number, the effective reproduction number takes the form of

If there is a single “average” infected individual at the onset of COVID-19 outbreak at time t = 0 then we have , as given in Eqs. (38), (42).

Parameters and initial conditions

For the parameters β1, β2, k and α that cannot be easily obtained, we use Markov Chain Monte Carlo (MCMC) simulations (Gamerman & Lopes, 2006) to estimate their values based on COVID-19 cases (https://covidmapper.ca/ge). The random-walk Metropolis-Hastings (MH) algorithm is used to adapt and readjust the MCMC procedure. We assume that the observations over time for data is Gaussian distributed (Gamerman & Lopes, 2006; Roda, 2020). For the parameter β3, it is about twice the parameter β2 (Li et al., 2020b; Moghadas et al., 2020). Note that transmission rates in this paper are real transmission probability at which a susceptible individual is infected by an infectious individual in a unit time. For the network topological parameters, the clustering coefficient φ and the initial value of average degree n0, we apply the method in Ref. (Luo et al., 2021). to calculate their values according to the distribution of household structure in the studied two cities (https://www.census.gov/da). The clustering coefficient φ is calculated by (2). The average degree is equal to the sum of the degrees of all nodes divided by the total number of nodes. At the initial time, the sum of the degrees of all nodes equals to , and the total number of nodes is . Thus, the initial value of average degree equals

In New York it appears that 88% of infected pregnant patients may be asymptomatic (https://www.washingtonpos, 2020). In Seattle more than half of the infected patients in a nursing home were asymptomatic (https://www.cidrap.umn.ed, 2020). Therefore, we consider that the probability of showing symptoms p equals to 0.25 for these two cities. For the remaining parameters, we obtain their values based on references (Luo et al., 2021; Moss et al., 2020; Park et al., 2020). All parameters are summarized in Table 2.

The initial values of variables in the pairwise model (3)–(35) are given as follows. The population sizes of New York and San Francisco are about 8, 336, 697 and 4, 264, 934, respectively (https://en.wikipedia.org/a; https://en.wikipedia.org/b; https://en.wikipedia.org/c; https://en.wikipedia.org/d; https://en.wikipedia.org/e; https://en.wikipedia.org/f). Let t0 stand for the starting date, which is March 27th, 2020. Values [E](t0), [A](t0), [SS](t0), [SE](t0), [SA](t0), and [SI](t0) are estimated by MCMC method. According to the collected data, the initial values of cumulative confirmed cases (denoted by C(t0)) and deaths caused by COVID-19 [R](t0) are obtained directly. The remaining initial values are obtained by data processing and calculations. Since the H individuals move into R compartment with proportion 1 − d or R compartment with proportion d, the initial value of recovered individuals with symptoms [R](t0) is calculated by the following approximation formula:

Since each E individual moves to I state with probability p or to A state with probability 1 − p, and p = 0.25, the ratio of the initial value of recovered individuals without symptoms [R](t0) and the initial value of recovered individuals with symptoms [R](t0) is less than (1 − p)/p = 3. We assume that the initial value of recovered individuals without symptoms [R](t0) is about 2.5 times of that of those with symptoms, that is,

For any time t, the number of cumulative confirmed cases is the sum of the confirmed symptomatic individuals [I](t), the number of recovered cases [R](t) and the number of died cases [R](t). Hence the initial value of the confirmed symptomatic individuals [I](t0) approximately equals to the difference between the number of cumulative confirmed cases and the sum of the number of recovered cases [R](t0) and the number of died cases [R](t0), that is,

Since the individuals in state I move into the H compartment with proportion ϵ, the number of the hospitalized individuals [H](t0) is calculated by the following approximation formula:

Quarantined individuals are caused by the diagnosis of infectious ones. The initial value of quarantined individuals approximately equals to the total number of close contacts of new confirmed cases in a incubation period: n0([I](t0 + 1/q) − [I](t0)). Among the number of quarantined individuals, [S0](t0) accounts for about 2/3, and the sum of [E0](t0), [A0](t0) and approximately accounts for 1/3 (https://www.soundhealthan). We assume that , then the approximation equations are

Since the number of all nodes is N which remains the same at any time, the initial value of unquarantined susceptible individuals [S](t0) is calculated by the following approximation formula:

For the initial values of links except estimated ones and [SR](t0), we determine their values based on the following approximation formula:where [X] or [Y] represents the number of nodes in state X or Y and [XY] represents the number of XY-links. According to formula (37), letting t = t0, we obtain the initial value of [SR](t0). All the initial values are summarized in Table 3.

Table 3

The initial values of variables.

Variable	New	Standard	San	Standard	Data
	York	deviation	Francisco	deviation	source
[S](t0)	8,158,859	–	4,259,357	–	(https://en.wikipedia.org/a; https://en.wikipedia.org/b; https://en.wikipedia.org/c; https://en.wikipedia.org/d; https://en.wikipedia.org/e; https://en.wikipedia.org/f)
[E](t0)	26,809	157.208 0	763.455 8	4.845 2	MCMC
[A](t0)	54,872	420.920 6	1,357	7.322 7	MCMC
[I](t0)	20,693	–	763	–	(46)
[H](t0)	4,139	–	152.7	–	(47)
[RA](t0)	11,285	–	431.7	–	(45)
[RI](t0)	4,514	–	172.7	–	(44)
[Rd](t0)	366	–	14	–	(https://covidmapper.ca/ge)
[S0](t0)	36,773	–	1,282	–	(48)
[E0](t0)	15,988	–	569.800 0	–	(49)
[A0](t0)	1,199	–	35.612 5	–	(50)
[R0A](t_0)	1,199	–	35.612 5	–	(51)
[SS](t0)	14,545,774	129,689	9,559,037	114,793	MCMC
[SE](t0)	66,189	266.492 9	1,764	27.537 7	MCMC
[SA](t0)	128,521	822.550 8	2,231	14.231 3	MCMC
[SI](t0)	2,063	21.625 0	1,102	12.209 7	MCMC
[SRA](t0)	13,910,916	–	5,207,559	–	(37)
[SRI](t0)	15,307	–	597.507 6	–	(53)
[EE](t0)	298.724 8	–	0.473 5	–	(53)
[EA](t0)	611.427 9	–	0.841 4	–	(53)
[EI](t0)	230.575 8	–	0.473 5	–	(53)
[ERA](t0)	125.745 3	–	0.267 7	–	(53)
[ERI](t0)	50.298 1	–	0.107 1	–	(53)
[AA](t0)	1,251	–	1.494 9	–	(53)
[AI](t0)	471.941 0	–	0.841 2	–	(53)
[ARA](t0)	257.374 7	–	0.475 7	–	(53)
[ARI](t0)	102.949 9	–	0.190 3	–	(53)
[II](t0)	177.973 8	–	0.473 4	–	(53)
[IRA](t0)	97.058 6	–	0.267 7	–	(53)
[IRI](t0)	38.823 5	–	0.107 1	–	(53)
[RARA](t0)	52.931 3	–	0.151 4	–	(53)
[RARI](t0)	21.172 5	–	0.060 6	–	(53)
[RIRI](t0)	8.469 0	–	0.024 2	–	(53)
n0	3.465 0	–	3.465 0	–	(43)

Results

Fitting results

Based on the diagnosed cases of COVID-19 infections (https://covidmapper.ca/ge), we adopt the MCMC method (Gamerman & Lopes, 2006) for 20 000 iterations with a burn-in of 10 000 iterations to fit the model (3–37) to estimate the parameters and the initial conditions of variables (see Table 2, Table 3). The transmission rates are subject to the condition that β1 ≤ β2. According to the shelter in place policy, t = 62 (May 28th) for New York, and t = 59 (May 25th) for San Francisco. Notice that there is no overfitting because the ratio of the number of data points (=115) and the number of estimated parameters or the initial values of variables (=10) is far over 5:1. The details have been presented in Section 2.4.

A good fitting between the model solution and real data in two cities is shown in Fig. 2. From the inset of Fig. 2(a), the real data point on April 23rd is a singularity with a negative value. This is due to a drop in the number of cumulative confirmed cases that day. In fact, ignoring this point has no effect on the overall fitting result. According to the estimated parameter values and initial conditions as given in Table 2, Table 3, we obtained the mean value of the reproduction number (95% CI of [1.444, 1.498]) in New York, (95% CI of [1.432, 1.459]) in San Francisco. Here CI indicates credible interval, the range containing a particular percentage of probable values. In other words, the 95% credible interval is simply the central portion of the posterior distribution that contains 95% of the values. In what follows, we mainly employ our model with estimated parameters to explore the epidemic behavior of COVID-19 infections in two cities.

Fig. 2

The fitting results of our model to real data of the cumulative and daily new cases (see inset) of COVID-19 infections in New York city and San Francisco city. The blue solid lines represent model simulation. Red dots and red lines are the real data. The grey area marks the 95% CI of MCMC estimations.

Fig. 3 shows the effective reproduction number. For these two cities, the effective reproduction number decreases quickly due to the social distancing. When shelter ended, the effective reproduction number increases quickly. For New York, the effective reproduction number is always smaller than one. But for San Francisco, the effective reproduction number is always larger than one, and the effective reproduction number is larger than 1.3 when shelter in place ended, which leads to a rapid increase of cumulative confirmed cases. Overall, the effective reproduction number is consistent with that shown on a popular website calculator for (https://rt.live/us/; https://rt.live/us/). The major difference is that when policy shifts the increase in our results is faster. The reason is that the intensity of intervention α is not equal to zero only at the intervention time. At the intervention time t, shelter ended, and the resumption of work and school began. This leads to an sudden increase in the total number of links in the network and hence an sudden increase in the average degree.

Fig. 3

The effective reproduction number.

Differences in COVID-19 infections

Perhaps surprisingly, under the assumption of no further changes in policy or transmission dynamics, the final size of the San Francisco epidemic is about 3.7 times that of New York, and that the epidemic duration of San Francisco is more than 140 days longer than that of New York. From Table 2, Table 3, transmission rates, population size and initial infections are different in New York and San Francisco. The transmission rate β1 is higher in New York than in San Francisco. The transmission rate β2 is slightly higher in San Francisco than in New York. According to Ref. (https://en.wikipedia.org/a; https://en.wikipedia.org/b; https://en.wikipedia.org/c; https://en.wikipedia.org/d; https://en.wikipedia.org/e; https://en.wikipedia.org/f), the population and area of New York are 8,336,697 and 785 square kilometers, respectively and those of San Francisco are 4,264,934 and 6405.5 square kilometers. So the population density of New York is about 16 times as dense as that of San Francisco. The starting date is not the true initial time of the disease, so it is of little significance to study the effect of the initial infections.

Here we only compare the effect of transmission rates and population density on disease dynamics by asking how the disease dynamics would progress if transmission rates or population densities were switched between the two cities. In detail, we set the transmission rates of New York being the estimated values of San Francisco to study the effect of transmission rates. According to the population sizes and areas in San Francisco and in New York, we calculate the population density of New York and modify the number of links [SS], [SE], [SA], [SR] accordingly. In a similar way, we investigate their effects in San Francisco. As shown in Fig. 4, population density has a major effect on disease transmission, and transmission rates have a smaller impact. For New York, when the population density drops to the same as San Francisco, the epidemic duration will be shortened by more than 140 days, and the final size will be reduced by more than 50%. When transmission rates in New York are made the same as those of San Francisco, the epidemic duration will be prolonged by more than 2 months, but the final size will be reduced by more than 11%. For San Francisco, when the population density increases to the same as New York, the epidemic duration will be extended by 19 days, and the final size will be increased 15-fold. When San Francisco transmission rates are made the same with those of New York, the epidemic duration will be shortened by 135 days, but the final size will be increased by more than 11%. Overall, the final size increases and the epidemic duration lengthens as population density increases. Also, as the rate of infection increases, so does the final size of the epidemic.

Fig. 4

The effect of transmission rates and population density on COVID-19 infections in New York city and San Francisco city. In panel (a), the green dot-dashed line is obtained by replacing transmission rates of New York with those of San Francisco, and the red dashed line is under the case of replacing population density of New York with that of San Francisco. In panel (b), the green dot-dashed line is obtained by replacing transmission rates of San Francisco with those of New York, and the red dashed line is under the case of replacing population density of San Francisco with that of New York.

Sensitivity analysis

To identify which parameters are important for the basic reproduction number R0, the final size and the effective reproduction number R(t), we conduct a sensitive analysis by evaluating the partial rank correction coefficients (PRCCs) for key parameters based on Latin hypercube sampling (Marino et al., 2008). All parameters are sampled 1000 times using a Latin hypercube sampling. For the probability density function (PDF) of these parameters, a normal distribution is selected for these parameters. Fig. 5 shows PRCCs which illustrate the dependence of R0, the final size and R(t) on all parameters, respectively. The panels (a) and (b) show PRCCs for R0. For New York, the four parameters β1, β2, q, φ have high correlations (|PRCCs| > 0.4) with R0. The most influential parameter on R0 is found to be the clustering coefficient φ (|PRCC| = 0.961 7), and then in descending order are transmission rate by A individuals β2 (|PRCC| = 0.882 3), transmission rate by E individuals β1 (|PRCC| = 0.735 6), incubation period for E individuals 1/q (|PRCC| = 0.660 5). The parameter γ2 has moderate correlations (0.2 ≤ |PRCCs| < 0.4) with R0. The other four parameters p, ε, γ1, γ4 have low correlations (0 ≤ |PRCCs| < 0.2) with R0. The PRCCs associated with these parameters except the parameter ε are statistically significant. For San Francisco, the five parameters β1, β2, q, φ, γ2 have high correlations with R0. The most influential parameter on R0 is found to be the clustering coefficient φ (|PRCC| = 0.982 2), and then in descending order are incubation period for E individuals 1/q (|PRCC| = 0.827 5), transmission rate by A individuals β2 (|PRCC| = 0.789 5), transmission rate by E individuals β1 (|PRCC| = 0.714 9), the recovery rate of A individuals γ2 (|PRCC| = 0.446 4). The other four parameters p, ε, γ1, γ4 have low correlations with R0. The PRCCs associated with these nine parameters are also statistically significant. From the panels (a) and (b), it is obvious to conclude that the higher the transmission rates are, the higher the basic reproduction number R0 is, and that the longer the incubation period is, the higher the basic reproduction number R0 is. In addition, there is a highly negative correlation between the clustering coefficient φ and the basic reproduction number R0. The shelter in place policy makes the clustering coefficient increase and then reduces the basic reproduction number.

Fig. 5

Partial rank correlation coefficients (PRCCs) for R0, final size and R(t) in New York and San Francisco. The panels (a) and (b) show PRCCs for R0 and each parameter, (c) and (d) show PRCCs for final size and each parameter, while (e) and (f) show temporal variation of the sensitivity of the effective reproduction number R(t) to each parameter. All parameters are sampled 1000 times using a Latin hypercube sampling. For the probability density function (PDF) of these parameters, a normal distribution is selected for these parameters.

The panels (c) and (d) show PRCCs for the final size. For these two cities, The PRCCs associated with the two parameters γ1 and γ4 are not statistically significant, thus we do not plot their PRCCs. For New York, this sensitivity analysis shows the first five parameters with most impact on the final size are transmission rate by A individuals β2 (|PRCC| = 0.615 8), the clustering coefficient φ (|PRCC| = 0.617 2), transmission rate by E individuals β1 (|PRCC| = 0.531 9), incubation period for E individuals 1/q (|PRCC| = 0.483 6), and probability of showing symptoms p (|PRCC| = 0.381 7). The PRCCs associated with the seven parameters (β2, φ, β1, q, p, γ2, ε) are statistically significant. For San Francisco, the first five parameters with most impact on the final size are probability of showing symptoms p (|PRCC| = 0.441 9), transmission rate by E individuals β1 (|PRCC| = 0.263 9), incubation period for E individuals 1/q (|PRCC| = 0.153 4), proportion of the confirmed symptomatic cases transferring to hospitalized cases ε (|PRCC| = 0.141 5), and transmission rate by A individuals β2 (|PRCC| = 0.098 1). The PRCCs associated with these five parameters are statistically significant. The PRCCs associated with the other two parameters φ and γ2 are not statistically significant.

The panels (e) and (f) show the temporal variation of the sensitivity of the effective reproduction number R(t) to each parameter. For these two cities, it is shown that variations in the clustering coefficient φ, the transmission rate by A individuals β2, transmission rate by E individuals β1, probability of showing symptoms p, and incubation period for E individuals 1/q dominate the PRCCs during the shelter in place policy. When the shelter in place policy ends, the correlations between these nine parameters (β2, β1, q, p, ε, φ, γ1, γ2, γ4) and the effective regeneration number R(t) decrease. During the shelter in place policy, the clustering coefficient φ is the most influential parameter on the effective reproduction number. There is a highly negative correlation between the clustering coefficient φ and R(t). The shelter in place policy makes the clustering coefficient increase and then reduces the effective regeneration number.

The effect of social distancing

Fig. 6 illustrates the effect of social distancing in New York and San Francisco. The panels (a) and (b) show the effect of social distancing on cumulative confirmed cases, while (c) and (d) show the effect of social distancing on the effective reproduction number. In all panels, black lines denote the actual case; green lines correspond to the case of keeping the shelter in place policy; red lines represent the case of changing the intensity of intervention α when shelter in place policy ended; blue lines represent bringing the shelter in place end time forward by one week; and purple lines correspond to delaying the shelter in place end time by two weeks. For New York, not lifting the shelter in place policy has a little effect on final size, but the epidemic duration will be shortened by more than 80 days. When the New York α keeps the same with that of San Francisco, that is, nearly doubling the intensity of intervention α, the final size will increase nearly 35%, and the epidemic duration will extend by one and a half more years. Bringing the New York shelter end time forward by one week will increase the final size by nearly 33%, and extend the epidemic duration by more than one year, whereas delaying the end time by two weeks has little effect on final size, but delays the epidemic duration by about 40 days. For San Francisco, not lifting the shelter policy will reduce the final size by 97%. When the San Francisco α keeps the same with that of New York, that is, reducing the intensity of intervention α by about half, the epidemic duration will be delayed by about 300 days and the final size will be reduced by about 21%; moving the shelter end time forward or later will shorten the epidemic duration and increase the final size. This means that New York will be less affected by reopening, but San Francisco will be more affected. Moreover, keeping a social distance helps reduce the final size of these two cities. The effective reproduction number in San Francisco is always more than one, which may lead to the results of San Francisco different from those of New York. It is noticed that these results are obtained under the assumption of no further changes in policy or transmission dynamics.

Fig. 6

The effect of social distancing in New York and San Francisco. The panels (a) and (b) show the effect on cumulative confirmed cases, while (c) and (d) show the effect on the effective reproduction number.

Discussion

Our comparative mathematical study of COVID-19 outbreaks in New York and San Francisco have allowed us to assess the impacts of population density and social distancing in these cities. Our approach was to use a pairwise network model to investigate COVID-19 epidemic in these two cities under a series of control measures. Using MCMC fitting to confirmed cases we were able to estimate the parameters and initial values of variables.

We explored the reasons for major differences between New York and San Francisco. Comparing the transmission in New York and San Francisco, we found that, under the current dynamics, the final size of San Francisco is predicted to be about 3.7 times as much as that of New York, and that the epidemic duration of San Francisco is more than 140 days longer than that of New York. This is because the effective reproduction number for New York was initially very high, but declined quickly due to social distancing. By way of contrast, San Francisco started with a lower effective reproduction number, but this did not decline as quickly. Unless there are major changes, the outcome for San Francisco may have an extended epidemic that ends off being larger than New York's.

We conducted a sensitive analysis by evaluating the PRCCs, and found that there is a highly negative correlation between the clustering coefficient and the basic reproduction number and the effective reproduction number. The shelter in place policy makes the clustering coefficient increase and then reduces the basic reproduction number and the effective reproduction number.

By mathematically switching population densities between San Francisco and New York, we were able to assess the impact of high-density interactions on COVID-19 dynamics. As expected, we found that population density has a major effect on disease transmission. The greater the population density is, the greater the final size of the epidemic is. Alternatively, by switching the transmission rates between the two cities we found that the final size increases as the rate of infection increases. As expected, social distance measures have a crucial and highly sensitive role to play in determining disease outcomes. Small changes such as increasing or decreasing the average number of contacts have the effect of allowing or preventing the outbreak in both New York and San Francisco. Thus, keeping strict social distance is conductive to disease prevention and control, and when social distance measures are relaxed, the epidemic may become more serious. Importantly, the timing of social distancing was also key, which is consistent with the conclusions in (Morris et al., 2020) on robust epidemic control.

Other approaches have been used to study the impact of physical distancing. In British Columbia, Canada a Bayesian model was used to analyze the population according to those willing and able to participate in distancing measures. The population was divided into two types of individuals: the non-physical-distancing and the physical-distancing (Anderson et al., 2020). By way of contrast, our model is based on the requirement of “shelter in place”. The quarantined individuals have close contacts with the confirmed individuals but do not participate in the COVID-19 infection.

As with many modelling problems, it is not always the case that a more complicated model is better. For example, using the Akaike Information Criterion (AIC) for model selection, Roda et al. showed that a susceptible-infectious-removed (SIR) model performs much better than a more realistic SEIR model in representing the information contained in the confirmed-case data (Roda et al., 2020). Of course, with the AIC model selection approach, the appropriate choice of model complexity depends on the type and amount of data available, and could change as more data become available. In the context of the models in this paper, assessing the appropriate level of complexity for pairwise interaction models, given the available data, is a possible venue for future work in this field, although it is likely that the optimal complexity from an information theoretic perspective will change as more data become available.

Funding

This work is partially supported by the grants 61 873 154 and 12 101 573, grants 2020XM18, Shanxi Provincial Department of Science and Technology COVID-19 Emergency Special Fund grants 202003D31011/GZ, and Fundamental Research Program of Shanxi Province grants 20 210 302 124 608 and 20 210 302 124 381, and partially supported by a (MAL), Discovery Grants (HW and MAL), Discovery Accelerator Supplement Award (HW), and an grant 202 100 502.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.