Optimizing agent-based transmission models for infectious diseases

Model structure and implementation

We have opted for a model structure consisting of households, schools, workplaces
and districts similar to published studies 6],10]. Figure 1 shows a schematic overview of the locations, which represent a group of people we
define as a “cluster”. Social contacts can only be made within a cluster. During nighttime,
people can have social contacts with members of their household and home district.
During daytime, people stay at home or go to a workplace or school depending on their
age, which also determines their day district. Contact between infectious and susceptible
people may lead to disease transmission, which is a stochastic process based on social
contact rates, infectiousness and susceptibility.

Fig. 1. Social contact structure. People are member of a household cluster and the corresponding
home district at night. During daytime, people can stay at home or go to a school
or workplace in a day district.

Figure 2 presents the model implementation with a general class diagram. We use a Simulator to organize the activities from the people in the Area. The Area has a Population, different Cluster objects and a Contact Handler. The Contact Handler performs Bernoulli trials based on the age of the contacts and random numbers. We
included a 2 ×2 social contact matrix, based on literature 20]-22], in which the transmission rate is doubled for contacts between children (18 y).
Each Cluster contains links to its members. The Population stores all person data (id, age, household, home district, day cluster, day district
and health related parameters) within or without Person objects but we elaborate further on this issue in the next paragraph. An infection
is assumed to follow a temporal pattern of Susceptible-Exposed-Infectious-Recovered
(SEIR) states similar to an influenza-like disease 6],10]. After infection, people need 2 days of latency (infected but not infectious) before
becoming infectious and 6 days to recover and acquire immunity against future infections.

Fig. 2. Model design: classes and compositions. The digits represent the number of links that
are possible. E.g., the Area can have 1 or many (*) Cluster objects, but a Cluster can only be part of 1 Area. The models differ in the implementation of the grey classes: FLUTE has less Cluster types in Area and the Population in SID does not contain Person objects.

We have constructed three implementations for the previous described transmission
model: “FLUTE” and “FRED” are based on the corresponding open source models and “SID”
has a novel data layout. The Area in FLUTE contains only home and day district Clusters. Membership to smaller sub-clusters like households, schools and workplaces can be
retrieved from stored cluster IDs in Person. People in a district that are also member of the same sub-cluster have two opportunities
for social contact and transmission. Therefore, during the processing of social contacts
in the district, sub-cluster IDs need to be checked. If two people from a district
are also member of the same sub-cluster, we used an aggregated transmission probability
instead of performing two random draws. The Area of FRED and SID has also separate households and day Clusters (= workplaces and schools). We illustrate the difference with the following pseudo-code
for the social contacts during nighttime with age dependent transmission probability
P _tr and P _tr? for one or two social contacts respectively:

Population structure

Data for an individual is stored as a Person object in FLUTE and FRED and the Population is a container of Person entities, stored consecutively in memory. In SID, the Population has a different container for each person attribute and the data of one person is
always located at the same index in each of those different containers. For example,
to access the age of person i in FLUTE or FRED we use “population[i].agepopulation[i].age” while in SID we use “population.ages[i]population.ages[i]”.

Population data have been extracted from the 2010 U.S. Synthetic Population Database
(Version 1) from RTI International 23],24] for Brooklyn and Nassau County, New York. Every county or state from this database
can be used to obtain individual age, household, school and workplace data. People
of 16 to 18 years of age with a school and work ID in the original database were assigned
to the school to guarantee that people were assigned to only one day cluster. To compare
different model implementations, we needed an extra social contact layer (Fig. 1). We have created home districts by adding households, sorted on ID, until a number
of 2000 people was reached. We assumed that household IDs are based on geographic
proximity and the threshold was adopted from Chao et al.6]. The day districts have been created analogously. The Nassau population consists
of 1.31 million people in 448 519 households and 140 861 day clusters. Brooklyn has
2.46 million people and the cluster sizes range from one up to 62 962 people. More
details on the study populations are listed in Table 1.

Table 1. Population statistics. Legend: [min – max] and (median)

The population data file determined the initial ordering of the person data in the
Population object. We used seven different orderings for the same population details: the original
sequence from the RTI database, a fully randomized order and population data sorted
according to household, day cluster, and both household (first) and day cluster (second),
and vice versa. To minimize the effect of random draws, we created five different
files for each ordering with a random component.

Algorithmic extension: sorting

The open source models 6],10] process disease transmission by looping over all members of a cluster and if a member
is infectious, to match them with all susceptibles. To reduce the total number of
operations, we introduced a modified algorithm in which the members of a cluster are
first sorted according health status before the infectious members are matched with
the susceptible members. A newly infected member is moved ahead of the first susceptible.
The member list obtains the following structure: First, recovered and infected (exposed
and infectious) members and second, susceptible members. The following pseudo-code
shows the sort algorithm for FRED and SID (the algorithm for FLUTE is structured analogously).

Parallelization: scheduling

The OpenMP API is often used for shared memory parallel programming in C/C++ 25]. In this programming model, subsets of a process are managed independently (=threads)
and share a global address space of a single or multiple processors which they read
and write asynchronously. For each cluster type (household, day district,…) in an
area, a person is a member of only one cluster. Therefore, clusters are stored per
type so that these containers can be processed in parallel without synchronization.
Parallel processing within one cluster would lead to synchronization overhead. The
workload distribution over the threads can be static or dynamic 25]. With static scheduling, a fixed number of tasks are assigned to each thread. In
dynamic scheduling, the workload is distributed over the idle threads until all tasks
are done. We have used workloads in chunks of one and ten clusters.

Inputs and work environment

We used a 2 × 2 transmission matrix and assumed that the transmission probability
(P _tr) is doubled for contacts between children (18y) 20]-22]. Similar to the literature 6],10], we estimated the relationship between P _tr and the basic reproduction number R ₀ by counting the number of secondary cases of one infected in a complete susceptible
population with seven P _tr values. Based on 4000 realizations with seven P _tr, we approximated R ₀ by exp(5507*P _tr-0.1911). The total run time depends on the clinical attack rate (AR, total fraction of the
population initially at risk that got infected) and for this reason, we performed
we performed benchmarks for a range of R ₀ values (1.1, 1.25, 1.4, 1.8 and 3). Each simulation was performed for 100 days. To
start the epidemic, we infected a random fraction of the population. After testing
seeding rates of 1e ^?2, 1e ^?3, 1e ^?4 and 1e ^?5, we observed limited impact on the number of cases for these ranges and selected
1e ^?4 as baseline setting.

We included the pseudo-random number generator (PRNG) from an open source software
package called TRNG 19],26], a portable and highly optimized library of parallelizable generators. To prevent
synchronization and latency, independent streams of random numbers are required for
each thread. We used the robust and versatile “leapfrog” method where the PRNG sequence
is distributed over p processes by calculating for draw i the i*(p-1)th number in the sequence. There are no recommendations to select PRNG seeds to obtain
different stochastic results, except that those seeds have to be different. Therefore,
the run index has been used to seed the PRNG.

An extended class diagram and the free open source code can be found in Additional
file 1 and Additional file 2 respectively. Additional file 3 contains a user manual to make use of the project software. During development, we
used the Google C++ Testing Framework 27] to perform detailed tests. These tests were applied in automated fashion with every
change in the code base via a continuous integration server 28]. The Templatized C++ Command Line Parser library 29] was used to transfer configurations to the executable. The project-code is standard
C++11 throughout, independent of external libraries and portable over all platforms
that have a GNU compiler (version 4.8 or later) available.

Timings presented in this paper were obtained from benchmarks on a cluster with Intel®;
Xeon®; E5-2680 v2 2.80 GHz CPU’s (release Q3’13) from the HPC core facility CalcUA
at the University of Antwerp. We confirmed our results with benchmarks on quad-core
Intel®; Xeon®; W5580 3.2 GHz (release Q1’09) CPU’s and AMD Opteron®; 6274 CPU’s. The
GNU compiler (4.8) was used in release mode with compiler optimization “-O3”. Additional
file 4 contains more info on the hardware and extra results. The open source tool PerfExpert
30] was used for profiling, as installed on the CalcUA cluster.

We performed additional benchmarks to explore the effect of cluster size, dynamic
clusters and increased model complexity on model performance. Methods and results
can be found in Additional file 5.