Matrix-Based Formulation of Heterogeneous Individual-Based Models of Infectious Diseases: Using SARS Epidemic as a Case Study.

Heterogeneities of individual attributes and behaviors play an important role in the complex process of epidemic spreading. Compared to differential equation-based system dynamical models of infectious disease transmission, individual-based epidemic models exhibit the advantage of providing a more detailed description of realities to capture heterogeneities across a population. However, the higher granularity and resolution of individual-based epidemic models comes with the cost of increased computational complexities, which result in difficulty in formulating individual-based epidemic models with mathematics. Furthermore, it requires great effort to understand and reproduce existing individual-based epidemic models presented by previous researchers. We proposed a mathematical formulation of heterogeneous individual-based epidemic models using matrices. Matrices and vectors were applied to represent individual attributes and behaviors. We derived analytical results from the matrix-based formulations of individual epidemic models, and then designed algorithms to force the computation of matrix-based individual epidemic models. Finally, we used a SARS epidemic control as a case study to verify the matrix-based formulation of heterogeneous individual-based epidemic models.

1. Introduction

Heterogeneity among the units composing a system is a very generic feature, which means that one never finds two units behaving in exactly the same way in most real systems [1]. In human societies, individual variability widely exists in social settings, cultural backgrounds, ages, genders, and occupations. Human beings have diverse features, physical characteristics, activities, and contact patterns. In epidemic outbreaks, individuals exhibit various health states, infectiousness, epidemic progresses, and behavior response. Despite these general facts, it is usually assumed that a population is homogeneous so as to build analytically tractable epidemic models, such as differential equation-based system dynamical models [2,3].

However, recent studies have shown that individual heterogeneity plays an important role in the epidemic spreading process [4,5,6]. The impact that individual heterogeneity has on epidemic diffusion has increasingly received more attention. Some researchers have concentrated on the effect of heterogeneous populations on epidemic dynamics [6,7,8]. They incorporated the heterogeneities of host populations as state variables, including age, stage, size [8], infection threshold [7], and resistant hosts [6]. Some other researchers have focused on epidemic models with spatial heterogeneities [9,10,11]. They utilized a patch structure where an epidemic system was divided into a number of patches to represent different spatial characteristics. Many researchers also investigated the impact of individuals’ heterogeneous behaviors and activities on epidemic diffusion [12,13,14,15]. Yang et al. [12] studied the impact of heterogeneous human activities on epidemic spread and found that the heterogeneity of activities affects spreading velocity remarkably. Merler et al. also studied how human motilities and population heterogeneity affect the course of an epidemic [13]. They found that the cumulative attack rate is positively correlated with average household size and the fraction of students in a population, and this was negatively correlated with the fraction of inactive population. Shang addressed the impact of awareness on epidemic outbreaks by proposing a mean-filed approach accommodating heterogeneous transmission rates [14]. He found that both contact awareness and local awareness can raise the epidemic threshold, while global awareness cannot. Cui et al. applied an edge-based compartmental model to study epidemic spreading dynamics with heterogeneous contacts [15]. In their models, contacts were grouped into two classes: strong contacts and normal ones. Additionally, the super spreading events that emerged in the severe acute respiratory syndrome (SARS) epidemic outbreaks of 2003 highlighted the necessity of a better understanding of social heterogeneities and individual-based models [16]. When studying super spreading events of the SARS epidemic using differential equation models, Mkhatshwa et al. [17] divided infectious individuals into two subgroups: super spreaders and regular spreaders. They then assigned a higher infection rate to super spreaders. Duan et al. employed heterogeneous and stochastic agent-based models to analyze the characteristics of infectious diseases’ super spreaders [18]. Researchers also paid attention to epidemic dynamics in heterogeneous networks with different connectivities [19,20,21], and went forward and investigated epidemic spread in weighted networks representing the heterogeneity of individual interaction strengths [22,23,24,25].

Differential equation models and individual-based models are the two most widely used approaches of representing epidemic diffusion [26]. The former makes some reasonable assumptions and simplifications so as to represent epidemic spread at a macro level [2]. A population is assumed to be homogeneous, well-mixed, and aggregated into several compartments according to people’s health states. The transitions of population between compartments are described by using differential equations with variables, such as the infection rate, onset rate of symptoms, and recovery rate. These assumptions and simplifications endow differential equation models with the advantage of performing theoretical analysis of macroscopic regularities of epidemic spread, such as the epidemic threshold and the final epidemic size. However, the advantage also comes with the limitation in representing epidemic diffusion in detail [27]. Firstly, the assumption of homogeneous and well-mixed population results in difficulties in representing the variants of individual microscopic attributes and behaviors. Secondly, a small set of variables parameterized with average quantities and mean values are inadequate to capture a variety of factors associated with the epidemic spreading process.

The individual-based model is a promising modeling paradigm to represent individual heterogeneities [28,29]. This approach can provide a more detailed depiction of realities to capture heterogeneity across individuals and incorporate the stochastic nature of infectious disease transmission. However, the higher granularity and resolution provided by individual-based models comes with the cost of increased computational complexities, growing computational power, and the development of intelligent computational algorithms [30]. The complexity of heterogeneous individual-based models leads to difficulty in mathematical formulation. Sometimes, heterogeneous individual-based epidemic models do not provide an insight into technical details, but only some qualitative representations. Consequently, we need to make great efforts to understand and reproduce existing individual-based epidemic models.

Our purpose is to realize the mathematical formulation of heterogeneous individual-based epidemic models, as well as make individual-based epidemic models much easier to be understood and reproduced. We used matrices and vectors to represent individuals, activity locations, social organizations, and relationship networks. Matrix-based representation of individual-based epidemic models reserves the capability of describing the heterogeneity of individuals and epidemic spreading process [8]. Furthermore, all of factors or elements related to epidemic diffusion are described as specific variables in matrices or vectors. According to the specific variables, it is feasible to build a mathematical formulation of individual-based epidemic models and algorithms. In addition, we can derive analytical results from the matrix-based representation of individual-based epidemic models. We applied the SARS epidemic as a case study to verify the matrix-based formulation of heterogeneous individual-based epidemic models.

2. Heterogeneous Individual-Based Models of SARS Epidemic

We here presented an improved agent-based model of the SARS epidemic based on our previous works, found in reference [18], which were used as a case study of matrix-based formulations of heterogeneous individual-based epidemic models. The main improvement of the model here is to increase the presentation of individuals’ daily commute activities among spatial locations into our previous works so as to study the heterogeneous spatial transmission pattern of the SARS epidemic. Consequently, we designed individuals’ daily activities at spatial locations, and introduced a bipartite network to describe individuals’ commute behaviors. We assumed a closed human society to study the spread and control of the SARS epidemic. We applied agent-based models to represent human individuals, spatial activity locations, social organizations in the closed society, and their behaviors. In addition, we employed complex networks to described the multi-relationships among these social entities. We designed individuals’ daily activities, such as sleeping, working, shopping, and recreation, which take place at different spatial activity locations. Agents commute among spatial locations to execute daily activities and contact each other. The contacts between susceptible individuals and infectious individuals cause the transmission of the SARS pathogen. Social organizations take some non-pharmaceutical interventions to mitigate the transmission of SARS epidemic, such as contact tracing, spatial activity location closure, individual quarantine, which may change the states of individuals and spatial activity locations. We then computed the discrete events and updated the states of social entities in time sequence to simulate the spread and control of the SARS epidemic.

2.1. Individuals’ Daily Commute Patterns

We represented individuals’ daily activities as mobility events, and randomly scheduled these discrete events in queuing models. We assumed individuals executed mobility activity events in each day.

We used a bipartite network with time-varying weights to model individuals’ daily commute patterns. There are two types of vertices in the bipartite network, including individual vertices and spatial location vertices. Network edges only exist between individual vertices and spatial location vertices. Edge weights denote the propensity that individuals visit spatial locations.

The bipartite network was designed in term of individuals’ living regularities. On account of the power-law characteristics of human dynamics [31], the edge weights of bipartite network were generated by a truncated power-law random variable described as where is a truncated power-law random variable. is a uniform random number. is the exponent of the power-law distribution. The probability that an individual moves to a spatial location in a commute activity is defined as where is the probability that individual i move to spatial activity location j. is the edge weight between individual i and spatial location j at time t. is the number of spatial locations individual i could visit at time t. It means spatial locations closed in the SARS epidemic outbreak is not allowed to be visited by individual i at time t.

2.2. Individuals’ Structured Contact Patterns

We also represented individuals’ contact activities as contact events, and randomly scheduled contact events in queuing models. According to survey results of human daily contacts [32], we used a normal random variable to represent the number of individuals’ daily contacts, which is described as where is the random variable of individuals’ daily contact number. and are both uniform random numbers in the range [0, 1]. In addition, considering the statistical features of human behaviors, we designed the time duration of each contact activity as a power-law waiting time distribution described as where is a random variable denoting the time duration of individuals’ contacts. is a uniform random number. is the exponent of the power-law distribution.

We built a contact network in term of individuals’ social relationships to represent individual contact patterns. The edge weights of contact network were generated by the formulation described as where is the edge weight between individual i and individual j. is the proportional parameter. and are the node degrees of individual i and individual j, respectively. Once individual i initiates a contact activity, it selects a neighbor in its personal contact network as a contact object according to the following rule where is the probability that individual i selects individual j as a contact object. here is the edge weight between individual i and individual j in the contact network. here is the number of neighbors in personal contact networks of individual i, which could be contacted by individual i. It means individuals who are quarantined or are not at the same spatial activity location cannot be contacted by individual i.

2.3. Infection Probability of SARS

SARS pathogens are transmitted from infectious individuals to susceptible individuals in their face-to-face contacts. Infectious individuals shed SARS pathogens through aerosols and droplets to susceptible individuals. Once susceptible individuals inhale a sufficient dose of SARS pathogens, they are infected by the SARS pathogen. Consequently, the infection probability of the SARS epidemic in a single contact is related to contact intensity and duration, infectivity of infectious individuals, and immunity of susceptible individuals. Here, we only considerer contact duration and infectivity. We assumed the infection probability of SARS epidemic in a single contact is where is the time duration of a single contact between a susceptible individual and an infectious individual. is the infectivity of the infectious individual. The infectivity is related to pathogen load and shedding rates of infectious individuals, and evolves along with the infectious period. According to consensus documents [33], we used a triangle distribution to model the infectivity (Inf) of infectious individuals, which is described as where and are the time points of beginning and ending infectious period. is the time point when infectious individuals have the max infectivity, which is set as day. According to Reference [18], we fitted the infection probability (Inf-Pro) of SARS in case infectious individuals have the max infectivity as where is the time duration of a contact between an infectious individual and a susceptible individual. We then designed infection probability by the proportion of infectivity with the max infectivity as where is the time of infectious period.

2.4. Heterogeneous Epidemic Progress of SARS

We assumed each infected individual went through four health state period, including susceptible period, latent period, infectious period, and recovered period. A susceptible individual is infected by the SARS epidemic in a single contact with an infectious individual. Then, the health state of the susceptible individual become latent. Several days later, the individual goes over the latent period, and become in infectious state. The infectious individual exhibits the capability of transmitting SARS pathogens to other susceptible individuals. On account of the onset of symptoms, infectious individuals will be admitted as patient cases of the SARS epidemic in hospitals. Admitted cases receives treatment in hospitals, and are quarantined until they recover from infectious state. Recovered individuals could be discharged from hospitals as they will not be infected by SARS epidemic again.

The duration of health state periods, admission time, and discharge time are fitted as Gamma distributions, of which the means (standard deviation) were 6.37 (16.69), 23.5 (62.1), 4.85 (12.19), and 23.1 (62.1) days, respectively [34]. To build the heterogeneous presentation of SARS progress, we applied Gamma random variables to scale the SARS epidemic progress.

3. Matrix-Based Formulation of the Models

To make heterogeneous individual-based models of SARS epidemic above much easier to be understood and reproduced, we applied matrices and vectors to describe objects and entities in the closed society. More technical details related to models and algorithms are represented in the matrix-based formulation. We firstly defined the models’ framework as a group of sets formulated as {, , , }. is the set of matrices representing entities and their attributes, such as individuals, spatial locations, health states. is the set of constraint conditions for entity matrices set . is the set of analytical results derived from . is the set of algorithms advancing the computation of heterogeneous individual-based epidemic models formulated by matrices.

3.1. Matrix-Based Formulation of Entities

3.1.1. Formulation of Individuals and Their Attributes

We assumed a closed population consisted of () individuals (Ind) in epidemic system. Individuals were represented as where is individual matrix with n rows and 1 column. is individual i in the population. The constraint condition of individual matrix is defined as , which means the row of individual matrix equals to the number of individuals in the population. Similarly, we used matrices to formulate individuals’ attributes, such as health states, epidemic progresses, infectiousness, and immunity.

We formulated individuals’ health states (HS) as where is health state matrix with n rows and h columns. h is the number of individual health states in epidemic progresses. means individual is in health state j (x) at time t. Susceptible state (S), latent state (L), infectious state (I), and recovered state (R) are marked as , respectively. The constraint conditions of health state matrix include , denoting the label of health state j of individual at time t, meaning individual have to be in only one health state at time t, and denoting all individuals are in one of health states at time t.

Epidemic progresses describe the transitions of individuals among different health periods. Individuals go through different periods of health states, which occupy different time durations, such as latent period, infectious period, and recovered period. We described epidemic progresses (EP) based on heterogeneous time scales as where is epidemic progress matrix with n rows and m columns. is the time duration of individual in the period of health state j (x). We described infection time (IT), latent period (LP), infectious period (IP), admission time (AT), and discharge time (DT) of individual as , , , , and , respectively. The constraint conditions of epidemic progress matrix include , meaning the time duration of health state period comes from a certain function, which is a stochastic Gamma distribution or a constant value, and meaning in case individual is not infected, all time durations of health state period are 0.

We applied a vector to represent the infectivity (Inf) of infectious individuals as where is infectivity matrix with n rows and 1 column. is the infectivity of individual at time t. The constraint conditions of infectivity matrix include , meaning the infectivity of infectious individuals is a function defined in Equation (8), and denoting individuals have infectivity only in case they are in infectious state.

Recovered individuals exhibit diverse abilities to resist pathogen. We formulated the immunity (Imm) of individuals as where is immunity matrix with n rows and 1 column. is the immunity of individual at time t. The constraint conditions of immunity matrix include , meaning the immunity of individuals is a certain function, and denoting individuals exhibit immunity in case they are in recovered state.

Once individuals are admitted as patients in hospitals, they are quarantined and prohibited to move outside and contact with others. We formulated individuals’ quarantine (Qua) state as where is quarantine matrix with n rows. is label denoting whether individual is quarantined. The constraint conditions of quarantine matrix include and .

3.1.2. Formulation of Spatial Locations and Their Attributes

Human mobility patterns depend on the transmission routes of infectious diseases. We assumed that a closed environment consists of many spatial locations. The number of spatial locations (SL) was assumed as . The spatial location matrix is formulated as where is spatial location i in the scene. The constraint condition of spatial location matrix is , which means the number of rows of spatial location matrix equals to the number of spatial locations in the closed environment.

During epidemic outbreaks, non-pharmaceutical interventions may be employed to close some spatial locations, such as campus, restaurants, and cinemas. We formulated the closure (Clo) state of spatial locations as where is closure state matrix with l rows. The constraint conditions of closure state matrix include and .

3.1.3. Formulation of Social Organizations and Their Attributes

We assumed there were social organizations (Org) in the closed environment, and formulated social organizations as where is social organization matrix with g rows. is social organization i. The constraint condition of social organization matrix is , which means the rows of social organization matrix equals to the number of social organizations.

Considering social organizations own spatial locations, such as office buildings, we formulated the ownerships between social organizations and spatial locations (Org-SL) as where is called as organization-location relationship matrix with g rows and l columns. is the label denoting whether social organization is sited in spatial location . The constraint conditions of organization-location relationship matrix including , , and . Furthermore, we also formulated the relationship between social organizations and individuals (Org-Ind) as where is organization-individual relationship matrix with n rows and g columns. is the label denoting whether individual belongs to social organization . The constraint conditions of organization-individual relationship matrix include , , and .

3.1.4. Formulation of Weighted Bipartite Networks

We applied a weighted bipartite network to represent individuals’ commute patterns among spatial locations. There are two types of vertices in the bipartite network, which denote individuals and spatial locations, respectively. Network edges only exist between individuals and spatial locations. Edge weights represent the propensity that individuals visit spatial locations. We formulated the bipartite network (BN) as where is bipartite network matrix with n rows and l columns. is the label denoting whether individual is linked to spatial location . The constraint conditions of bipartite network matrix include , , and . The edge weights of bipartite network (BNW) are formulated as where is the edge weights matrix of bipartite network with n rows and l columns. is the edge weight between individual and spatial location at time t. The constraint conditions of edge weights matrix of bipartite network include , , meaning the edge weights of bipartite network are generated by a function defined in Equation (1), and meaning the edge weight is 0 in case individual is not connected to spatial location .

Individuals’ current spatial positions (CSL) at time t are recorded as where is individual current spatial location matrix with n rows and m columns. is the label denoting whether individual is at spatial location at time t. The constraint conditions of individual current spatial location matrix include , , , meaning individual have to stay at only one spatial location, and meaning all individuals stay at the spatial locations.

3.1.5. Formulation of Weighted Contact Networks

We represented individuals’ heterogeneous contact patterns by using a weighted contact network, where nodes stood for individuals, edges represented individuals’ relationships, and edge weights described individuals’ interaction strength. We formulated the contact network (CN) as where is contact network matrix with n rows and n columns. is the label denoting whether individual is connected to individual . The constraint conditions of contact network matrix include , , and . Edge weights of the contact network (CNW) were formulated as where is the edge weight matrix of contact network with n rows and n columns. is the edge weight between individual and individual at time t. The constraint conditions of edge weight matrix of contact network include , meaning the edge weights of contact network is generated by a certain function defined in Equation (5), and meaning the edge weight between two individuals is 0 in case they are not connected with each other. We marked who were contacting with whom as where is contact object (CO) matrix. is the label denoting whether individual is contacting with individual at time t. The constraint conditions of contact object matrix include , , and meaning individuals cannot contact with themselves.

3.1.6. Formulation of Transmission Networks

We recorded the transmission networks of epidemics to analyze transmission routes of epidemics and derive the reproduce number of infected cases. The epidemic transmission network (TN) is represented as where is the label denoting whether individual is infected by epidemic from individual . The constraint conditions of transmission network matrix include , , and meaning individuals cannot transmit epidemics to themselves.

3.1.7. Formulation of Individuals’ Mobility Activities

We represented the heterogeneous mobility patterns of individuals, and scheduled individuals’ mobility activities as discrete events in queue models. Individual mobility event (ME) matrix is formulated as where is mobility event matrix with n rows and k columns. is the mobility event j of individual , which denotes individual will move to spatial location at time . The constraint conditions of mobility event matrix include meaning mobility events are scheduled in time sequence, meaning the time point when mobility event (METP) happens is generated by a certain function, meaning which spatial location individual moves (MESL) to at time is generated by a certain function defined in Equation (2).

3.1.8. Formulation of Individuals’ Contact Activities

We also represented individual heterogeneous contact patterns and scheduled contact activities as discrete events in queue models. Individual contact activity matrix is formulated as where is contact event (CE) matrix with n rows and k columns. is the contact event j of individual , which denotes individual initiates a contact event with individual at time . Furthermore, the contact activity sustains a duration time . The constraint conditions of contact event matrix include meaning contact events are scheduled in time sequence, meaning the time point (CETP) when contact events happen is generated by a certain function, denoting individuals select contact objects (CEO) according to a function defined in Equation (6), and meaning the duration time of contact activities (CED) is generated by a certain function defined in Equation (4).

3.2. Statistical Properties of Epidemic Diffusion

To conduct theoretical analysis of epidemic spreading regularities, we derived analytical results from matrix-based formulation of entities above. The numbers of susceptible, latent, infectious, and recovered individuals at time t were formulated as

The number of new infections (NI) at time t is derived as

The accumulated number of infected individuals (AI) is derived as where is the accumulated number of infections at time t. Then, we can get final epidemic size as .

To analyze the spatial patterns of epidemic diffusion, we derived the distribution of infectious individuals at spatial locations (SL-I) as where is the proportion of infectious individuals at spatial location at time t. is the label denoting whether individual is in infectious state and stays at location at time t. is the total number of infectious individuals at time t. Furthermore, we derived the density of infectious individuals at spatial locations as where is the density of infectious individuals at spatial location . is the number of individuals who are at spatial location at time t. Similarly, we got the distribution of infectious individuals at social organizations (Org-I) as where is the proportion of infectious individuals at social organization at time t. is the label denoting whether individual is in infectious state and belongs to social organization . The density of infectious individuals at social locations is formulated as where is the density of infectious individuals at social locations . is the number of individuals who belong to social organization . We derived the reproduce number of epidemics according to transmission network matrix as where R is the reproduce number of epidemics. is the number of infected individuals who produce more than one of the second-generation cases.

3.3. Algorithms of Computation

We designed algorithms to force the computation of heterogeneous individual-based epidemic models formulated by matrices and vectors. The computational framework of our models is illustrated in Figure 1.

Figure 1

The computational framework of algorithms.

There are four computational processes invoked in time ticks, which contain eight major modules composed of algorithms in the framework. These modules are described as follows.

Computing engine: an engine that forces the computation of models. It advances the computation of models in time order. In each time tick, the engine invokes the computing of mobility activity progress, contact activity progress, and epidemic progress. Furthermore, it realizes time management service that calculates date list. Besides, it invokes events schedules to generate and update mobility events and contact events at the beginning of each day.

Mobility events schedule: a container that generates, updates, stores mobility events. At the beginning of each day, it generates new mobility events for each individual to update event schedule.

Contact events schedule: a container that generates, updates, store contact events. At the beginning of each day, it generates new contact events for each individual to update event schedule.

Mobility activity progress: a progress that executes the computation of a mobility activity including choosing expected location, calculating individual’s position, etc. A weighted spatial network is applied to build the model of individual mobility.

Contact activity progress: a progress that executes the computation of a contact activity including choosing expected contact object, calculating contact duration, calculating infectious probability, etc. A weighted contact network is applied to build the model of individual contact.

Epidemic progress: a progress that represents the evolution of health states of infected individuals. It fulfills the computation of individual transitions among different health states by using the infected time and the duration of each health state.

A weighted temporal-spatial network: a weighted network that combined a spatial network with a contact network. The weighted temporal-spatial network is applied to build the model of individual mobility and contact behaviors.

Random number generator: a generator that produces random numbers subject to various distributions. These random numbers are utilized to represent the stochastic nature of epidemic spreading and individual behaviors (e.g., infection probability, mobility time, contact time, contact duration, the durations of health states, etc.).

These modules implement different services in order to execute the computation of epidemic models. It is essential to design algorithms to realize the capabilities of providing these services in these modules. We summarize four major algorithms represented as follow.

3.3.1. Algorithm of Computing Engine

Computing engine implements three jobs described in Figure 2. Computing engine firstly initializes models and parameters described as matrices and constraints. Then, computing engine enters the main loop of calculations and advances the computation of models in time ticks. In each time tick, computing engine invokes epidemic progress, mobility activity progress, and contact activity progress. Meanwhile, computing engine record the state of epidemic diffusion and update mobility and contact event schedules. Computing engine finally outputs analytical results and terminates the computation of models.

Figure 2

Algorithm of computing engine.

3.3.2. Algorithm of Mobility Activity Progress

We represent mobility activity progress in Figure 3, which contains three jobs: querying a mobility activity in mobility event schedule; selecting an activity location according to weighted spatial mobility network; and updating the current location matrix.

Figure 3

Algorithm of computing mobility activities.

The function in Figure 3 decides which spatial location is selected as individuals’ destination in their daily mobility activities according to the probability formulated as where is the probability that individual selects spatial location as his destination in a daily mobility activity. The probability is defined by using the proportion of the edge weights of weighted bipartite network. denotes spatial location is not closed in epidemic outbreaks. This equation is a matrix-based formulation of Equation (2).

3.3.3. Algorithm of Contact Activity Progress

We described contact activity process in Figure 4, which contains four jobs: querying a contact activity in contact event schedule; selecting a contact object according to weighted contact network; updating current contact object matrix; and computing the infection at the end of contact activities.

Figure 4

Algorithm of computing contact activities.

The function in Figure 4 decides who are going to contact with whom in individuals’ contact activities according to the probability described as where is the probability that individual selects individual as a contact object when individual initiates a contact event. The probability is defined by using the proportion of edge weights of weighted contact network. denotes individual is not quarantined at time t. ensures individual and individual are simultaneously at the same spatial location at time t. This equation is a matrix-based formulation of Equation (6).

The function decides whether susceptible individual is infected by epidemics from infectious individual in a contact activity for a time duration . The infection function is formulated as where is the time point when the infectious period of individual begins. is the time point when the infectious period of individual ends. (days) is the time point when the infectivity of infections individual reaches the maximum. This equation is a matrix-based formulation of Equation (10).

3.3.4. Algorithm of Epidemic Progress

We represented epidemic progress in Figure 5, which contains three jobs: checking the health state of individuals according to health state matrix ; judging whether the current health state of individuals is over according to epidemic progress matrix ; and then updating health state matrix.

Figure 5

Algorithm of computing epidemic progress.

4. Simulating the Spread and Control of SARS Epidemic

To test and verify matrix-based formulation of heterogeneous individual-based transmission models of the SARS epidemic, we used the matrix-based models and algorithms to simulate the spread and control of the SARS epidemic.

4.1. Experiment Design

We assumed a closed society containing individuals and spatial activity locations. We used the method of randomly adding edges to build a bipartite network and a contact network. The edge weights of bipartite network were generated by a truncated power-law random variable . The exponent of power-law distribution was set as 1.5 in Equation (1) according to refence [31]. The edge weights of contact network were generated according to Equation (5) whose proportional parameter was set as 0.1, and scaling parameter was set as 1.

We assumed each individual executed five mobility activities in a day. Individuals’ daily contact number approximates to a normal distribution with mean 8.6 (standard deviation 4.25) [34]. The time duration of contact activities was represented as a truncated power-law random variable in the range minute with an exponent 1.5 [31].

We studied the effect of non-pharmaceutical interventions on the SARS epidemic diffusion, including contact tracing and workplace closure. In case susceptible individuals ever contacted with infectious individuals were found out, they were quarantined. Furthermore, once quarantined individuals had the onset of symptoms, they were admitted to be infected by the SARS epidemic. The admission time of traced infectious individuals was changed to be a Gamma distribution with a mean (standard deviation) of 1 (0.25) day. In case workplaces were closed, individuals stayed at home or had activities in other spatial locations. The changes were realized by revising bipartite network matrix, contact network matrix.

4.2. Experimental Results

We randomly chosen an individual as the infection source of the SARS epidemic, and executed the epidemic control policies of contact tracing and workplace closure at various time in simulation experiments. We conducted a set of 100 experiments for each group of control policies, and illustrated experimental results in Figure 6.

Figure 6

Experimental results. The different days in the legends of subfigures denote the days until workplace closure and contact tracing were implemented. (a) The time evolution of the density of infectious individuals. (b) Final epidemic sizes with different days until contact tracing and workplace closing policies were implemented. (c) The distribution of the density of homes via the number of infected individuals. (d) The distribution of the density of workplaces via the range of infected individuals.

We observed the evolution of the density of infectious individual over time in Figure 6a. The slower we executed a workplace closing policy, the higher peak the curves got. Moreover, the curves quickly drop down after the peak, and may have serious oscillations in the second half of time. In Figure 6b, we found that with the increase of time when we began to use workplace closing policies, the data points denoting the final epidemic size in experiments scattered at the positive direction of Y axis. The final epidemic size has a higher mean value. In Figure 6c,d, we found that the densities of home places and workplaces approximated to power distributions over the number of infected individuals and the range of infected individuals, respectively. Furthermore, in case we adopted a workplace closing policy earlier, the epidemic invaded a fewer families and workplaces.

We found the correlation effect caused by contact activity heterogeneity and activity location heterogeneity in Figure 6. Individuals’ contact events are randomly scheduled in a single day according to a uniform distribution. Moreover, individuals conduct five mobility events in a single day and have activities at different locations in the same moment. When executing a contact event, an individual is only able to contact a neighbor that is also in his current activity location. So, the heterogeneities in the beginning time of contact events, the time duration of contacts, the contact probabilities between neighbors, the beginning time of mobility events, and the probabilities of activity locations bring up temporal-spatial constraints for individuals’ contact behaviors. At the same time, the contact patterns of individuals lead to the temporal-spatial characteristics of epidemic diffusion. Due to workplace closure interventions, the epidemic diffusion across different workplaces is mitigated. It is shown in Figure 6d that workplace closure interventions stop most workplaces from epidemics. However, individuals are likely to stay at home when they cannot go to workplaces. The increase of contact behaviors taking place at home enhances epidemic spread among families. This explains why the tails of curves raise up in Figure 6c.

5. Conclusions

Individual heterogeneity is a key factor bridging the link between micro behaviors and macro phenomena. As a bottom-up method, individual-based epidemic models exhibit the capability of representing the heterogeneities in epidemic spreading process. However, the complexity of individual-based models results in the limitation of mathematical formulation.

In this paper, we used matrices to formulate heterogeneous individual-based epidemic models. We adopted matrices and vectors to represent entities and their attributes in epidemic systems, such as individuals, individuals’ health states, the infectivity of infectious individuals, spatial locations, and bipartite networks. Then, we derived analytical results of epidemic diffusion from the matrix-based formulation of entities, and designed algorithms to force the computation of heterogeneous individual-based epidemic models. Finally, we applied the SARS epidemic as a case study to test and verify the heterogeneous individual-based epidemic models in matrix-based formulation. In this case study, we investigated the effect of non-pharmaceutical interventions on the SARS epidemic control, including contact tracing and workplace closure. Experimental results indicate that in case of executing workplace closing interventions early, we can obtain a better effect on mitigating SARS epidemic diffusion. Moreover, we also analyzed the correlation effect caused by contact activity heterogeneity and mobility location heterogeneity, as well as the temporal-spatial patterns of epidemic spread.

Heterogeneous individual-based epidemic models formulated by matrices and vectors are scalable and flexible. We could add more matrices into the models to describe more entities in epidemic systems. Furthermore, we could improve the algorithms or design new algorithms to expand the capability of the models. More analytical results could be formulated from matrix-based representations of entities and their attributes. In addition, heterogeneous individual-based epidemic models formulated by matrices and vectors not only can capture individual heterogeneity, but also have the ability of performing theoretical analysis of epidemic diffusion.