Inferring HIV Transmission Network Determinants Using Agent-Based Models Calibrated to Multi-Data Sources

Abstract

(1) Background: Calibration of Simpact Cyan can help to improve estimates related to the transmission dynamics of the Human Immunodeficiency Virus (HIV). Age-mixing patterns in sexual partnerships, onward transmissions, and temporal trends of HIV incidence are determinants which can inform the design of efficient prevention, and linkage-to-care programs. Using an agent-based model (ABM) simulation tool, we investigated, through a simulation study, if estimates of these determinants can be obtained with high accuracy by combining summary features from different data sources. (2) Methods: With specific parameters, we generated the benchmark data, and calibrated the default model in three scenarios based on summary features for comparison. For calibration, we used Latin Hypercube Sampling approach to generate parameter values, and Approximation Bayesian Computation to choose the best fitting ones. In all calibration scenarios the mean square root error was used as a measure to depict the estimates accuracy. (3) Results: The accuracy measure showed relatively no difference between the three scenarios. Moreover, we found that in all scenarios, age and gender strata incidence trends were poorly estimated. (4) Conclusions: Using synthetic benchmarks, we showed that it is possible to infer HIV transmission dynamics using an ABM of HIV transmission. Our results suggest that any type of summary feature provides adequate information to estimate HIV transmission network determinants. However, it is advisable to check the level of accuracy of the estimates of interest using benchmark data.

1. Background

The Human Immunodeficiency Virus (HIV) infection remains a major public health concern in sub-Saharan Africa (SSA). The vast majority of new HIV infections are caused by unprotected sexual intercourse [1]. Thus, sexual transmissions are believed to be the main driver of the epidemic in that part of the world [2,3]. High-risk sexual behaviour of individuals, such as having a large number of concurrent and lifetime sexual partners, and a high frequency of unprotected sex increases the risk of acquiring HIV infection [2,4,5,6,7]. Thus, for HIV and other sexually transmitted diseases, individual behaviours affect the transmission dynamics [8]. Our understanding of patterns of age-related sexual partner choices at population level known as age-mixing patterns [9] can help to undertake behaviour change interventions to prevent HIV transmission.

Temporal trends of HIV incidence [10], help us to see the impact of interventions on reducing new infections among populations [11,12]. Given the type and quality of ongoing HIV interventions, new infections may (or may not) produce secondary HIV transmissions. Thus, temporal trends of incidence alone do not give a full picture of the impact of interventions. Given a time interval, in a population where one infected individual generate only one secondary infection, if many of them are infected in that population, the number of new infections may be the same as when few infected individuals generate many secondary infections. Hence, we need also to estimate the distribution of onward HIV transmission.

Therefore, from a public health point of view, estimating age-mixing patterns in sexual partnerships, the distribution of onward HIV transmission, and temporal trends of HIV incidence, are important to inform the design of more effective and efficient HIV prevention and linkage-to-care programs.

Most of traditional HIV studies look at prevalence, incidence, distribution of the infection in different populations, and the risk of acquiring the infection. These types of studies use mostly epidemiological data, clinical, and sexual behavioural survey data [2,13]. Either models are fitted to these data [14,15], or statistical analysis is performed on them [3,13,15,16]. However, if we want to understand some key features of HIV infection, such as mixing between different groups of the populations of interest, the rapidity of viral spread within populations, transmission clusters and where new infections occur, and the transmission of drug resistance [17,18], we concur with the idea of undertaking phylogenetic analysis which uses viral sequence data. As HIV is transmitted across a sexual network, the structure and dynamics of the sexual network shape the evolution of the virus population [19]. The analysis of HIV phylogenetic trees (constructed from HIV sequences sampled from infected individuals) has the ability to determine estimates of some transmission network determinants of realised sexual networks [20].

The main output of a phylogenetic analysis is a phylogenetic tree, which is a diagram that represents evolutionary relationships among organisms or genes [21]. If constructed from viral sequences of infected individuals, it holds information about the structure of a sexual network [19], and the transmission dynamics [22,23]. The shape of a phylogenetic tree consists of two properties, the topology (particular branching pattern for a tree), and the branch lengths [24]. Thus, some topological features from a phylogenetic tree, such as the Colless imbalance and “ladder length”, are properties that have been identified to distinguish different transmission patterns [25]. In a “super-spreader” pattern, most infected people produce no or one onward transmission, while a small fraction produce a high number of onward transmissions [25]. In chain-like transmission, the vast majority of infected people produce exactly one onward transmission, while a few exceptions produce two or three onward transmissions. In a homogeneous transmission pattern, onward transmissions are Poisson-distributed. Thus, summary features from phylogenetic tree data were used to calibrate modelling framework and estimate HIV transmission networks [26], and sexual contact networks [19]. Despite some challenges related to phylogenetic analysis assumptions [17], viral sequence data were used to infer temporal trends of HIV incidence [22,27,28], as well as the age difference in HIV transmission clusters [18], and other important epidemiological parameters such as the basic reproductive number $R_0$ [29].

The estimates of the abovementioned HIV transmission network determinants (age-mixing patterns in sexual partnership, HIV temporal trend of incidence, and onward transmission) give an overview of the extent of age-mixing patterns in the sexual network, HIV transmission dynamics, and the effect of interventions. These estimates may be elusive or difficult to estimate, and we often rely on models which we fit to existing data. Fitting models to empirical data is known as calibration, and it entails searching the best fit values for the models’ parameters which can generate the data, given the observed summary features (known as target summary statistics) from observed data. We used Simpact Cyan [30], an integrated modelling framework which is an agent-based model (ABM) used to run a simulation study from which we explore its calibration with regard to estimation of determinants of the HIV transmission network. Thus, since this study is a simulation study, the observed data would be the benchmark data generated with a known vector of parameter values. This has been elaborated in detail in the methodology section.

Some authors argued that combining different data sources such as epidemiological data, social network data, phylogenetic tree data, and many more, has the potential to improve epidemic estimates, and thereafter advance public health interventions [31,32,33]. Therefore, we hypothesised that combining summary features from epidemiological and behavioural data, and phylogenetic tree data to calibrate the Simpact Cyan model [30] can improve estimates for age-mixing patterns in sexual partnerships, the distribution of onward transmissions, and the temporal trends of HIV incidence.

Simpact Cyan, being an HIV simulation tool, can therefore be applied and utilised as a modelling framework which simulates sexual dynamic networks, HIV transmission, and viral evolution across transmission networks. It has potential to help us understand some parts of complexity of sexual partnerships across sexual networks. Furthermore, it helps to advance methodological studies related to the use of HIV viral data that aim to understand HIV transmission networks [34]. This simulation study opens research pathways for calibration studies which would use Simpact Cyan, focusing on HIV inferences by shedding light on key considerations such as those from phylogenetic trees for summary features and parameter estimation. Thus, with this simulation study, we wanted to see if, by using synthetic benchmarks data, we can infer HIV transmission network determinants using an agent-based model simulation tool for HIV infection transmission.

2. Methods

2.1. Simulation Tool: Simpact Cyan

In this simulation study, we used Simpact Cyan [30], a modelling framework which simulates demographic events (birth and deaths), sexual dynamic networks, HIV transmission across sexual networks, and viral evolution across transmission networks. The core of the simulation tool is based on agent-based models, and the interactions of agents for sexual partnerships and HIV transmission, together with other events, are defined by hazard functions. These hazard functions are mathematical quantities which give the rate at which an event such as establishing or dissolving a sexual partnership, HIV transmission, birth, or death occurs. With fine-grained detail of agent-based models, we were able to add another component to simulate viral evolutionary dynamics above the transmission network. This allowed for the use of multiple data sources including sequence data in our simulation and analysis. More details of the components of the simulation tool are described in the Appendix A. The details of the Simpact Cyan simulation model are also given in the Appendix A, which is a part of the paper’s supporting information of the paper [34].

2.2. Simulation of HIV Epidemic

We simulated sexual partnership dynamics, HIV transmission, and viral evolution across transmission networks using Simpact Cyan [30,35]. The simulation time was 40 years (0 to 40 years). Due to the phylogenetic component of the simulation (the age of the root sequence we used), we considered calendar time, and the simulation time became 40 years between 1977 and 2017. We defined an initial population of men (5000) and women (5000). During the simulation, a sexual partnership was set to start 10 years earlier before the introduction of HIV in the population to allow equilibrium of partnership establishment/dissolution events. To initiate HIV transmission, 10 individuals aged randomly between 20 and 50 years were seeded with being HIV positive allowing the HIV model to run its course. Thus, HIV infection was introduced in 1987, and subsequent transmissions occurred up to 2017.

The simulation framework contains different hazard functions that capture the dynamic of sexual partnerships, HIV transmission, and ART intervention, among others. These functions were established based on different assumptions to ensure the simulated dynamics of HIV epidemic can be a proxy, hence depicting the real world situation. For instance, HIV transmissions are assumed to occur as a function of partner HIV status, viral load levels, antiretroviral therapy (ART) intervention and follow-up of the HIV positive partner. The gradual treatment eligibility based on the Cluster of Differentiation 4 (CD4) protein counts as implemented in the real world [36,37], was considered in the epidemic simulation. We considered 100, 150, 200, 350, 500, 700 CD4 counts for 23, 25, 28, 33, 36, 39 simulation years, respectively.

At the 40 years of simulation time point (2017 in calendar time), the simulated HIV epidemic exhibited different trends of prevalence between women and men. Younger women had higher prevalence compared to men of the same age groups as we can see in Figure 1. Between 39 and 40 years of simulation time (2016 and 2017 in calendar time), younger women (below 25 years) had higher incidence compared to men of the same age group as we can see in Figure 2. Thus, the simulated epidemic mimics commonly observed trends in generalised HIV epidemic settings in SSA [38,39,40].

The viral evolution across the transmission network was simulated using a transmission tree built from the transmission network, to simulate viral sequences (each per individual) with a general time-reversible $\left(GTR\right)+\Gamma$ model in Seq-Gen with a seed sequence [41]. The main advantage of this integrated simulation framework is that it enables estimation/computation of known measurements associated to the HIV epidemic, which may be elusive or hard to estimate in real life.

Key assumptions for HIV simulation and molecular evolution of the simulation model were provided. For HIV simulation, we had considerations on demographics, sexual partnerships, HIV transmission, ART interventions, and disease progression. For molecular evolution of the virus, we considered population-level variation of viral sequences with GTR substitutions of nucleotides. More details on these assumptions mentioned above are given in the Appendix A (

Table A2. Main assumptions for the simulation of the dynamic sexual partnership network and HIV transmission.

Sub-Component	Assumptions
Demographic	Birth: when there is a sexual partnership formation, a conception event will be scheduled; after a conception event is triggered, a new birth event will be scheduled, so that the woman in the relationship will give birth to a new person at a specific time, and the gender will be determined by the boy/girl ratio. Mortality: normal mortality model follows Weibull distribution, and the time for AIDS mortality was determined as the time of infection plus the survival time.
Sexual partnership	We considered sexual partnerships such that the preferred age gap differed from one person to the next, but there was also an age dependent component in this preferred age gap, and we allowed for the weight of the age gap terms to be age-dependent. We assumed the age gap to be normally distributed. The hazard function for partnership depended on the number of partners the man and woman in the relationship had. The debut age was set to be 15 years for men and women. Once a sexual partnership was established, it was subject to dissolution as well.
HIV transmission	Transmission is likely to occur when one individual in a sexual partnership is infected with HIV. Transmission depends on the viral load level of an infected individual. A woman’s susceptibility to HIV infection was considered to depend on her age. Set-point viral load is the viral load that the person has during the chronic stage. In the acute stage or in the Acquired immunodeficiency syndrome (AIDS) stages, the configuration values $person.vsp.toacute.x$, $person.vsp.toaids.x$ and $person.vsp.tofinalaids.x$ cause the real viral load to differ from the set-point viral load in such a way that the transmission probability is altered: the hazard for transmission will increase by the factor x that is defined in this way. We also assumed that, once an infected individual is on ART, s/he cannot transmit the infection; thus, we set $monitoring.fraction.lo{g}_{v}iralload$ to 0.
ART interventions	If this CD4 count is below the threshold set in $monitoring.cd4.threshold$, the person will be offered antiretroviral treatment (ART). Depending on the person’s willingness to accept treatment, treatment will then be started. We considered gradually increasing ART eligibility based on CD4 count thresholds as implemented in real life.
Disease progression	Follow up on the progress of the disease is performed by inspecting the person’s CD4 count. When a person receives treatment, the viral load is lowered and if the person drops out of treatment the viral load will increase.

and

Table A3. Parameters for the HIV 1 substitution $GTR+\Gamma +I$ model.

Name	Value
Relative Frequencies
Adenine (A)	0.3857
Cytosine (C)	0.1609
Guanine (G)	0.2234
Thymine (T)	0.2300
Relative substitution rates
r(A → G) = r(G → A)	2.9114
r(A → C) = r(C → A)	12.5112
r(A → T) = r(T → A)	1.2569
r(G → C) = r(C → G)	0.8559
r(G → T) = r(T → G)	12.9379
r(C → T) = r(T → C)	1.0000
Rate heterogeneity
Shape parameter	0.9
Number of gamma rate categories	4
Fraction of invariant sites
Proportion of invariant sites (I)	0.5230
Evolutionary rate
Substitutions/site/year 2	$4.75\times {10}^{-3}$

² which is the branch length scaling factor.

), which describes in more detail the simulation model and its components.

To generate benchmark data of reference, we assigned values to 13 selected parameters which are some of the known parameters to have direct effects on sexual partnerships, and HIV transmission; more details for those parameters can be seen in the Appendix B. With those parameter values, we ran the integrated simulation framework 1120 times in order to take into account stochasticity. The simulation framework with values of the selected parameters was considered a “master model”. From that model, we computed the measurements of interest in our study which were parameters for age-mixing patterns in sexual partnerships, temporal trends of HIV incidence, and the distribution of onward transmissions. We also computed two types of summary features: epidemiological and behavioural summary features, and phylogenetic tree summary features. These summary features were used to calibrate the default integrated simulation framework in three scenarios: (i) with epidemiological and behavioural features, (ii) with phylogenetic tree features, and (iii) with combined features.

In our simulation, for a man i with n partnerships, there were n age gap preferences, this gave clustered data regarding age gaps. The equation below (1) explains the variation of the man’s age gap preference for his women partners with a linear mixed random-effect model [42]:

$$y_{ij}={\beta }_0{x}_{ij}+{\beta }_1+b_i+{ϵ}_{ij}$$

(1)

where $y_{ij}$ represents the age gap preference of woman j in sexual partnership with man i, and $x_{ij}$ was the predictor which was the age of the man i. The parameters ${\beta }_0$ and ${\beta }_1$ represent the fixed effects, while $b_i$ parameters represent the random effects.

In our calibration exercise, we had two types of parameters to estimate age-mixing patterns in sexual partnerships; one was a set of descriptive parameters for age difference between a male and his female partners, namely the average age difference (AAD) across sexual partnerships, and the standard deviation of the age difference (SDAD). We also had parameters from the linear mixed random-effect model (LMM) given by Equation (1), which was fitted to the clustering data of age gap preference for men. Fitting the LMM to the data using the R package $lme4$ [43], we obtained other age-mixing parameters: the within-subject standard deviation of age differences (WSD), the average variation of age gap within the clusters of men’s age gaps; the between-subject standard deviation (BSD), the average variation of age gap between the clusters of men’s age gaps; the slope and intercept which show the overall population level trend of age gap preference.

For distribution of onward transmissions, we computed the mean, median, and standard deviation of the number of infected individuals, secondary generated by an infected individual who died between 35 and 40 years simulation time.

In addition for temporal trends (35–40 years of simulation time) of HIV incidence, we computed incidence values per year per age group (15–24, 25–29, 30–34, 35–39, 40–44 and 45–49 years) among men and women. We considered one year interval for incidence between 35 and 40 years of simulation time.

The calibration process entails the use of summary features (from benchmark data) which are also known as target statistics—which are always considered to be measurable from real data, and calibrate the default simulation framework (Simpact Cyan). Given the benchmark summary features (target statistics) and a parameter space where we can search parameter values, the calibration process aimed to find best guess of parameter values which produce summary features which are close to those from benchmark data. We then use these parameter values to run the simulation framework, and compute the parameters for age-mixing patterns in sexual partnerships, temporal trends in HIV incidence, and the distribution of onward transmissions at a certain level of accuracy when we can compare them to what we obtained in the benchmark data.

All the above mentioned transmission network determinants and summary features were computed from data sampled between 35 and 40 years of simulation time, which was 2012 and 2017 in calendar time.

2.3. Selected Parameters for Calibration

Simpact Cyan has many parameters, as we can see in its documentation [35]. Details about all events, associated settings, and their parameters are explained in the documentation. In our calibration exercise, we selected parameters of sexual partnership events with related settings, sexual partnership dissolution events, and HIV transmission events. More details of these parameters can be found in the Appendix A, where we described the Simpact Cyan simulation model, in the Appendix A.2 for events’ hazard functions and associated settings.

2.4. Selected Summary Features

The summary features were chosen with regard to how we believe they can capture information of the transmission network determinants (age-mixing patterns in transmission, onward transmission, and temporal trends of incidence). The benchmark values for epidemiological and behavioural features were computed from demographic, sexual behavioural, HIV epidemic, and intervention benchmark data. The phylogenetic tree features were computed from a phylogenetic tree built from sequence data. The first calibration process was conducted by considering only features from epi-behavioural data; the second was conducted with phylogenetic features; and the third was conducted with combined features. More details about summary features are given in the Appendix B.

The following summary features were selected for epidemiological and sexual behaviour data: (i) population growth rate between 35 and 40 years of simulation time; (ii) HIV infection point prevalence (40 years of simulation time) for men and women in different age groups (15–24 years, 25–29 years, 30–34 years, 35–39 years, 40–44 years, and 45–49 years) at 40 years of simulation time (end of simulation); (iii) HIV infection incidence (39–40 years of simulation time) for men and women in different age groups (15–24 years, 25–29 years, 30–34 years, 35–39 years, 40–44 years, and 45–49 years); (iv) prevalence of concurrent partnerships (39.5–40 years of simulation time); (v) relationship per person per year (35–40 years of simulation time); (vi) mean, median, and standard deviation of age difference between infected individuals (35–40 years of simulation time); (vii) anti-retroviral treatment coverage (from 33 to 40 years of simulation time); and (viii) viral load suppression at 40 years of simulation time for individuals with ≥15 years but <50 years.

On the other hand, we derived the following summary features from phylogenetic tree data: (i) mean, median, and standard deviation of of internal nodes; (ii) mean, median, and standard deviation of tree height; (iii) Colless and Sackin indexes; (iv) minimum, maximum, mean, median, standard deviation, 1st and 3rd quantiles of branch length; (v) mean, median, standard deviation of tips depths; (vi) mean, median, standard deviation of of nodes depth; and (vii) maximum height. Whenever phylogenetic data were used, the calibration process was conducted with summary features computed at 100% of sampling coverage in a scenario where we assumed that missing sequences were missing completely at random (MCAR).

Since, the master model was executed 1120 times, to summarise the summary features to be used during calibration, we considered mean values for each target statistic. We did the same regarding mean values for parameters of the HIV transmission network determinants.

2.5. Calibration Scheme

The calibration process was performed using Latin Hypercube Sampling (LHS) [44,45], a technique used to generate parameter(s) values in given parameter(s) space(s), together with the rejection method of Approximation Bayesian Computation (ABC) [46], an approach used to choose the best candidate of parameter(s) values given the summary features from benchmark data. Because we had to retain only 30% best fit parameter values with ABC calibration, we had to be able to run the model after calibration with 1120 parameter combinations; thus, we had to generate 3752 parameter values. With the Latin Hypercube sampling technique, we set 3752 sampling cells for parameter space for each parameter among the 13 parameters. Each parameter had 3752 values, which made a $3752\times 13$ matrix of parameter combinations. After running the default simulation framework 3752 times, we retained 30% of best parameter values using Approximation Bayesian Computation on each of the 13 parameters. Thus, each parameter had a vector of parameter values.

In order to conduct the calibration in a uniformly defined parameter spaces, we considered a closed set (which is a metric space) around each parameter $p_i$, with radius $p_i/2$, which is characterised by low values $(p_i-p_i/2)$ and high values $(p_i+p_i/2$) if $p_i$ is positive or vice versa if $p_i$ is negative. For all 13 parameters with known values of the benchmark scenario, we set 13 parameter spaces, $\Omega$ (vectors of parameter spaces), one per each parameter, for calibration:

$$\Omega =\left\{\begin{array}{cc}{\Omega }_{min}=P-\frac{1}{2}P\hfill & \hfill and{\Omega }_{max}=P+\frac{1}{2}P:P\ge 0\\ {\Omega }_{min}=P+\frac{1}{2}P\hfill & \hfill and{\Omega }_{max}=P-\frac{1}{2}P:P<0\end{array}\right.$$

(2)

where P is the vector of values of the 13 parameters used to generate the benchmark data. The calibration process was performed thrice, since we had three scenarios defined by summary features. With these retained parameters, we ran the simulation framework and computed the transmission network determinants.

We explored two options to run the simulation framework after calibration: the first one was to consider the vectors of selected parameter(s) and make several parameter combinations and run the integrated simulation framework. For example, if each of the 13 parameters had 400 parameter values, we had 400 parameter combinations, and we ran the simulation framework 400 times. The second option was to take the median of the parameter values for each of the 13 parameters, and with this unique parameter combination run the simulation framework 1120 times. Thus, for each calibration scenario, we computed the transmission network determinants twice. After that, we had to measure the goodness-of-fit (GoF) by comparing obtained values to what we obtained in the benchmark data. We used mean relative error (MRE) to evaluate the GoF. Thus, in this study, the two main components of model-calibration were: the parameter-search strategy by LHS and the ABC, and the GoF measure by the mean relative error [14].

The work-flow of this calibration simulation study was subdivided into three main steps, as we can see in Figure 3:

Step 1: with specific parameter values, we generated a benchmark data set, from which we computed the age-mixing in sexual partnerships, onward HIV transmission, and temporal trends of incidence, and values of summary features from epidemiological, behavioural, and phylogenetic tree data.

Step 2: we reset the modelling framework at default parameter values, and with the given summary features from the benchmark data, coupled with our prior knowledge about plausible parameter space, we tried to produce the best fit parameter values which could produce summary features close to those from the benchmark data.

Step 3: we took new parameter values from calibration as input to the modelling framework and ran the modelling framework to compute age-mixing in sexual partnerships, onward HIV transmission, and temporal trends of incidence.

In terms of computation time, the first and third steps take relatively the same amount of computational time, whereas the second step takes much more time due to the parameter search for the best fit.

The performance of the calibration method, and the modelling framework, was evaluated by computing the relative error between outputs from the benchmark and after calibration, including benchmark parameters and those from calibration, and associated summary features. In the following Results section, we discuss the accuracy of age-mixing in sexual partnerships, onward HIV transmission, and temporal trends of incidence, after the three calibration scenarios. Details on accuracy of retained parameters during calibration scenarios, and associated summary features, are given in the Appendix B.

All scripts to reproduce the results and data generated are publicly available at a GitHub repository (https://github.com/niyukuri/calibration_sympact, accessed on 15 November 2020).

3. Results

The following results are the mean values for age-mixing patterns in sexual partnerships, temporal trends in HIV incidence, and the distribution of onward transmissions. As stated earlier, after each calibration scenario, we had two options to run the simulation framework, either using all retained parameter values or their medians; we present results from both options.

3.1. Age Mixing Patterns in Sexual Partnerships

Compared to the benchmark values, the post-calibration estimates of the descriptive statistics AAD, and SDAD, were better estimated in all scenarios, as we can see in

Table 1. Mean values for estimates of age-mixing patterns in sexual partnerships from the master model and after calibration with features from HIV epidemiology and sexual behaviour (Calibration 1), phylogenetic tree (Calibration 2), and a combination of the two types of features (Calibration 3).

Estimate	Benchmark	Calibration 1	Calibration 1 *	Calibration 2	Calibration 2 *	Calibration 3	Calibration 3 *
AAD	13.782	13.861	14.167	13.151	13.111	13.633	13.83
SDAD	6.15	6.114	6.17	6.008	6.113	6.073	6.164
BSD	2.138	2.013	2.077	2.026	2.137	2.004	2.098
WSD	1.812	1.781	1.774	1.837	1.807	1.811	1.779
Slope	0.315	0.297	0.302	0.316	0.319	0.305	0.305
Intercept	−2.378	−2.111	−2.336	−2.099	−2.2	−2.164	−2.264

${}^{*}$ Values computed when we run the simulation with medians of the retained parameter values during calibration.

. The same is true for BSD, WSD, the slope, and intercept of the linear mixed-effects model.

Mean values without * were computed with vectors of the retained parameter values during calibration, and those with * were computed with medians of the retained parameter values. This is also applicable for the results of other HIV transmission network determinants in subsequent subsections.

From the mean relative error in

Table 2. Mean relative error for the average of age mixing estimates compared between benchmark data and after calibration with features from HIV epidemiology and sexual behaviour (MRE 1), a phylogenetic tree (MRE 2), and a combination of the two types of features (MRE 3).

Estimate	MRE 1	MRE 1 *	MRE 2	MRE 2 *	MRE 3	MRE 3 *
AAD	0.16399	0.03948	0.19238	0.05177	0.16766	0.0325
SDAD	0.10248	0.04145	0.10447	0.04005	0.10367	0.04208
BSD	0.17373	0.06907	0.16063	0.06254	0.16258	0.06986
WSD	0.09965	0.057	0.10925	0.05368	0.10193	0.05984
Slope	0.20558	0.08395	0.21915	0.07315	0.20482	0.07606
Intercept	0.37034	0.16361	0.39138	0.16511	0.39999	0.16004

${}^{*}$ Error obtained when we run the simulation with medians of the retained parameter values during calibration.

, overall, a two-by-two comparison of outputs obtained with vectors of retained parameters (calibration columns without *), and those with median values of retained parameters (calibration columns with *) for each calibration scenario, the second option performed well. Thus, after calibration, the median values of retained parameter values gave the best results of estimates of age-mixing patterns in sexual partnerships.

Values of mean relative error without * were computed by comparing benchmark outputs and those after calibration with vectors of the retained parameter values, and those with * were computed by comparing benchmark outputs and those after calibration with medians of the retained parameter values. This is also applicable for the mean relative error of other transmission network determinants.

As we can see in the error Table 2, a two-by-two difference check among relative errors (e.g., MRE 1 and MRE 1*) for different scenarios showed that the error varies between tenths and hundredths. Focusing on outputs produced with the parameters’ median values (columns 3, 5, and 7 in Table 2), which seems to have less error values. The comparison of the results after the three calibration scenarios showed that the difference between benchmark values and those obtained after calibration was around a few hundredths, except for the intercept which was around few tenths. When looking at parameters one by one, overall, the decreasing trend of accuracy was observed as follows: AAD, SDAD, WSD, BSD, Slope, and Intercept. Thus, the slope and intercept of the LMM were estimated with relatively low accuracy compared to other parameters.

3.2. Onward Hiv Transmission

For onward HIV transmission distribution,

Table 3. Mean values for estimates of HIV onward transmissions from the master model and after calibration with features from HIV epidemiology and sexual behaviour (Calibration 1), a phylogenetic tree (Calibration 2), and a combination of the two types of features (Calibration 3).

Estimate	Benchmark	Calibration 1	Calibration 1 *	Calibration 2	Calibration 2 *	Calibration 3	Calibration 3 *
Mean	2.242	2.384	2.245	2.323	2.24	2.356	2.246
Median	1.391	1.351	1.31	1.391	1.303	1.354	1.337
Standard deviation	2.12	2.622	2.27	2.428	2.296	2.579	2.273

${}^{*}$ Values computed when we run the simulation with medians of the retained parameter values during calibration.

shows that, compared to benchmark values, the mean, median, and standard deviation of the number of onward transmissions were better estimated in all scenarios. It also shows that, on average, one infected individual would produce more than two other infections or less, given the standard deviation value.

A two-by-two comparison of outputs obtained with vectors of the retained parameters (calibration columns without *), and those with median values of the retained parameters (calibration columns with *) for each calibration scenario of the mean relative error in

Table 4. Mean relative error of average values of onward HIV transmission statistics compared between benchmark data and after calibration with HIV epidemiology and sexual behaviour (MRE 1), a phylogenetic tree (MRE 2), and a combination of the two types of features (MRE 3).

Estimate	MRE 1	MRE 1 *	MRE 2	MRE 2 *	MRE 3	MRE 3 *
Mean	0.17513	0.1349	0.15716	0.12686	0.16807	0.13272
Median	0.3271	0.31999	0.34116	0.31962	0.32561	0.32279
Standard deviation	0.46932	0.32561	0.39235	0.32676	0.45539	0.33333

${}^{*}$ Error obtained when we run the simulation with medians of the retained parameter values during calibration.

showed that the second option performed well. Thus, after calibration, the median values of retained parameter values gave the best results of distribution of onward HIV transmissions.

In the error Table 4, a two-by-two difference check among relative errors (e.g., MRE 1 and MRE 1*) for different scenarios showed that the error varies around tenths. However, the outputs produced with the parameters’ median values (columns 3, 5, and 7 in Table 4) had less error values. When looking at parameters one by one, overall, the decreasing trend of accuracy was observed as the following: mean, median, and standard deviation.

3.3. Temporal Trend of Hiv Incidence

The benchmark trend of incidence in the five year period (2012–2017) for men and women in different age groups (in the first column of Figure 4 or Figure 5) showed that it was decreasing in all gender and age groups. Younger women between 15 and 24 years had the highest incidence compared to other age groups in women and men. The decreasing trend meant that ART interventions had effects on HIV transmission.

After the three calibration scenarios, the incidence data points computed with vectors of retained parameters were presented in Figure 4. The comparison between the three scenarios and the benchmark (columns of Figure 4) showed that the temporal trend of incidence was over-estimated in almost all calibration scenarios. However, much more over-estimation was after calibration of the simulation framework with phylogenetic tree features, followed by the scenario of combined features. The same conclusion was drawn for Figure 5, where incidence values were computed after running the simulation framework with median values of retained parameter values after calibration.

As previously seen for other estimates, a two-by-two comparison of incidence obtained with vectors of the retained parameters (calibration columns without * in

Table 5. Mean relative error for average incidence values in male populations compared between benchmark data and after calibration with HIV epidemiology and sexual behaviour (MRE 1), phylogenetic tree (MRE 2), and a combination of the two types of features (MRE 3).

Year	Age Group	MRE 1	MRE 1 *	MRE 2	MRE 2 *	MRE 3	MRE 3 *
2013	15–24	1.05475	0.88048	1.34137	1.0745	1.06295	0.93416
	25–29	0.87811	0.76712	0.95202	0.88461	0.89938	0.78916
	30–34	0.91606	0.8875	0.96849	0.91955	0.96005	0.88494
	35–39	0.78408	0.66974	0.7441	0.67547	0.75629	0.67397
	40–44	0.89173	0.78456	0.79473	0.73747	0.91572	0.74142
	45–49	1.23131	1.04045	1.13016	0.99237	1.22711	1.00805
2014	15–24	1.02293	0.89883	1.26296	1.0851	1.08533	0.97876
	25–29	0.96249	0.88694	1.10231	1.02698	0.98951	0.90741
	30–34	0.86408	0.87018	0.90373	0.89681	0.92165	0.84121
	35–39	0.97827	0.87733	0.93042	0.88425	0.99837	0.8737
	40–44	1.06991	0.91003	0.92059	0.91656	0.97718	0.88609
	45–49	1.45055	1.29848	1.27103	1.16657	1.4448	1.15002
2015	15–24	1.15855	0.89525	1.33417	1.25347	1.27041	1.03238
	25–29	1.00525	0.93552	1.09395	1.16055	1.12176	0.9917
	30–34	0.92342	0.89038	0.91207	0.96075	0.94996	0.93358
	35–39	1.07437	0.93358	0.95423	0.90022	0.90746	0.87915
	40–44	1.07316	0.88081	0.96276	0.80236	1.00017	0.85782
	45–49	1.75576	1.50713	1.53888	1.44127	1.55175	1.37258
2016	15–24	1.13425	0.98986	1.33453	1.22738	1.25674	1.0544
	25–29	1.02155	0.96523	1.23536	1.05684	1.05984	0.96505
	30–34	0.91918	0.94391	0.97993	1.03122	0.97103	0.97724
	35–39	1.18662	1.2673	1.25378	1.12695	1.16183	1.20835
	40–44	1.03665	0.96006	0.96548	0.83396	0.90925	0.92091
	45–49	1.47005	1.41002	1.40545	1.29973	1.38828	1.31173
2017	15–24	1.35979	1.1209	1.35127	1.29633	1.39386	1.21792
	25–29	0.95806	1.01259	1.08813	1.01433	0.99	0.97441
	30–34	1.00497	0.97803	1.08916	1.01595	1.02405	1.00899
	35–39	1.34227	1.46893	1.3267	1.30232	1.35038	1.37524
	40–44	1.13972	1.07222	1.20909	0.99736	1.05262	1.03022
	45–49	1.58792	1.51293	1.74517	1.38875	1.47887	1.41031

${}^{*}$ Error obtained when we run the simulation with medians of the retained parameter values during calibration.

and

Table 6. Mean relative error for average incidence values in the female population compared between benchmark data and after calibration with HIV epidemiology and sexual behaviour (MRE 1), a phylogenetic tree (MRE 2), and a combination of the two types of features (MRE 3).

Year	Age Group	MRE 1	MRE 1 *	MRE 2	MRE 2 *	MRE 3	MRE 3 *
2013	15–24	0.53162	0.47872	0.52974	0.49637	0.49674	0.46769
	25–29	1.1711	1.10191	1.02598	1.08174	1.06856	1.07796
	30-34	1.4554	1.59738	1.48093	1.4674	1.57672	1.38211
	35–39	1.27389	1.56912	1.43407	1.45138	1.4723	1.39072
	40–44	1.5237	1.77582	1.58568	1.65677	1.69882	1.47681
	45–49	1.41214	1.95796	1.45311	1.63971	1.51047	1.46651
2014	15–24	0.54265	0.50475	0.58463	0.5347	0.54235	0.48523
	25–29	1.35538	1.30714	1.25242	1.19749	1.26711	1.21816
	30–34	1.58765	1.74194	1.51293	1.54455	1.59564	1.3997
	35–39	1.2958	1.44317	1.48171	1.56704	1.47045	1.27727
	40–44	1.47626	1.65227	1.57011	1.61447	1.45799	1.44798
	45–49	1.51765	1.7599	1.51771	1.5789	1.49514	1.36147
2015	15–24	0.57955	0.52997	0.60027	0.57768	0.58239	0.53849
	25–29	1.51821	1.57128	1.53055	1.49222	1.59013	1.42023
	30–34	1.46685	1.64833	1.52767	1.55935	1.57546	1.48934
	35–39	1.87122	1.85779	1.71533	1.56098	1.75373	1.34148
	40–44	1.52558	1.72044	1.53849	1.80755	1.46566	1.42572
	45–49	1.72515	1.60911	1.48315	1.52458	1.62834	1.40109
2016	15–24	0.58743	0.60129	0.61815	0.62563	0.59805	0.60696
	25–29	1.46132	1.5059	1.41493	1.4317	1.48655	1.40186
	30–34	1.52139	1.58994	1.67728	1.57346	1.53513	1.44627
	35–39	1.67303	1.88848	1.78009	1.90138	1.91832	1.4736
	40–44	1.40186	1.73526	1.42221	1.50055	1.78263	1.34671
	45–49	1.52093	2.17204	2.05912	1.89981	1.45132	1.52321
2017	15–24	0.70051	0.6883	0.79729	0.74652	0.75446	0.71317
	25–29	1.58232	1.67152	1.58225	1.48288	1.50227	1.45087
	30–34	1.87383	2.07842	1.77046	1.9108	1.87073	1.54384
	35–39	1.73521	1.61932	1.6716	2.07031	1.49161	1.6285
	40–44	1.94166	2.25749	1.89354	1.8269	1.85843	1.5933
	45–49	1.26963	1.61916	1.24791	2.0685	1.32776	1.45465

${}^{*}$ Error obtained when we run the simulation with medians of the retained parameter values during calibration.

) and with median values of the retained parameters (calibration columns with * in Table 5 and Table 6) for each calibration scenario; we can see that the second option performed well although it does not give the best results of distribution of the temporal trend of HIV incidence.

If we consider mean relative error after calibration when the simulation was performed by median values of retained parameters (columns 4, 6, and 8) for both Table 5 and Table 6, the comparison between the error tables showed that men had higher error values compared to women. For the female population, younger women (15–24 years) had the lowest error values (Table 6).

When we compare the temporal trend of incidence between men and women (rows of Figure 4) for the benchmark and after calibration (in all three scenarios), we can see that incidence values for both populations had the same descending trend with time. The same conclusion was drawn from Figure 5, where incidence values were computed after running the simulation framework with median values of retained parameter values. However, when looking at the relative error, the latter did not decrease or increase with time. Thus, the over-estimation of incidence values did not vary with time, but across different age groups and gender as we can see in Table 5 and Table 6.

4. Conclusions

The objective of this simulation study was to investigate if accuracy of estimates obtained by calibration of the Simpact Cyan simulation framework with combined summary features can be improved. We expected that the estimates of HIV transmission network determinants (age-mixing patterns, distribution of onward HIV transmission, and temporal trend of HIV incidence) could be improved by fitting the simulation framework to combined summary features. From the values of mean relative errors in the three calibration scenarios, we have seen that we produced relatively the same outputs after calibration of the simulation framework by epidemiological and behavioural features, phylogenetic tree summary features, and combined features. Thus, any scenario can help to calibrate the simulation framework. These outputs results were supported by the parameter values retained in calibration for all scenarios (

Table A4. Median values and inter-quantile range for retained parameter values during calibration with epidemiological and behavioural features (Calibration 1), phylogenetic tree features (Calibration 2), and combined features (Calibration 3); together with input values used to run the master model (benchmark data), and the search parameter space for calibration.

Parameter	Calibration 1	Calibration 2	Calibration 3	Input	Space
$dissolution.alpha.0$	−0.514 [−0.654, −0.397]	−0.52 [−0.652, −0.39]	−0.514 [−0.653, −0.391]	−0.52	[−0.78, −0.26]
$dissolution.alpha.4$	−0.053 [−0.064, −0.041]	−0.052 [−0.063, −0.04]	−0.053 [−0.064, −0.041]	−0.05	[−0.075, −0.025]
$formation.hazard.agegapry.baseline$	1.957 [1.513, 2.401]	1.976 [1.546, 2.409]	1.965 [1.54, 2.398]	2	[1, 3]
$person.agegap.man.dist.normal.mu$	10.277 [8.287, 12.324]	9.666 [7.433, 11.913]	10.11 [7.949, 12.195]	10	[5, 15]
$formation.hazard.agegapry.gap.agescale.man$	4.898 [3.84, 6.128]	4.989 [3.728, 6.248]	4.96 [3.839, 6.191]	5	[2.5, 7.5]
$person.agegap.man.dist.normal.sigma$	0.255 [0.2, 0.308]	0.244 [0.186, 0.299]	0.252 [0.198, 0.306]	0.25	[0.125, 0.375]
$formation.hazard.agegapry.numrel.man$	−0.299 [−0.371, −0.23]	−0.302 [−0.37, −0.237]	−0.302 [−0.371, −0.239]	−0.3	[−0.45, −0.15]
$formation.hazard.agegapry.numrel.diff$	−0.101 [−0.126, −0.076]	−0.099 [−0.122, −0.074]	−0.1 [−0.124, −0.075]	−0.1	[−0.15, −0.05]
$hivtransmission.param.a$	−1.016 [−1.215, −0.739]	−0.977 [−1.195, −0.726]	−0.997 [−1.216, −0.739]	−1	[−1.5, −0.5]
$hivtransmission.param.b$	−85.888 [−109.426, −65.1]	−84.888 [−107.133, −64.703]	−85.927 [−109.915, −65.038]	−90	[−135, −45]
$hivtransmission.param.c$	0.589 [0.501, 0.662]	0.599 [0.512, 0.676]	0.59 [0.507, 0.662]	0.5	[0.25, 0.75]
$hivtransmission.param.f1$	0.048 [0.036, 0.061]	0.048 [0.036, 0.06]	0.048 [0.036, 0.06]	0.048	[0.024, 0.072]
$hivtransmission.param.f2$	−0.142 [−0.177, −0.106]	−0.142 [−0.18, −0.106]	−0.141 [−0.177, −0.105]	−0.14	[−0.21, −0.07]

). When we look at these as shown in the Appendix B, we observe that all retained parameter values were closer to those used to generate benchmark data; the error between them was at the order of ten thousandths (Figure A1).

If we multiply error values by 100, for all calibration scenarios, AAD, SDAD, BSD, and WSD estimates of age-mixing patterns have less than 7% of mean relative error. However, the slope and intercept of the mixed effect model for age-mixing patterns had 8 and 16% of mean relative error, respectively. For the distribution of onward HIV transmission, the error for the mean was around 12%, whereas for the median and standard deviation it was around 30%. For the temporal incidence trend, the majority of incidence values had more than 100% of relative error, which made incidence the poorly depicted measurement after calibration.

Therefore, estimates of age-mixing patterns in sexual partnerships had relatively lower relative error (around 10%) when compared with the benchmark values and those obtained after calibration in any calibration scenario. They were followed by the distribution of onward HIV transmission.

Studies on sexual partnerships and age difference between partners showed that age gap in sexual partnerships and subsequent HIV transmission could be the main driver of the HIV transmission cycle between generations in some settings [18,47]. One of the studies used sero-conversion data from a population-based cohort of men and women [47], and another one used phylogenetic pairings from transmission clusters computed from a phylogenetic tree [18], both in KwaZulu Natal. Thus, as we have shown, Simpact Cyan can be calibrated with epi-behavioural and phylogenetic tree data, providing a technical platform to understand the effect of age-mixing patterns on HIV transmission dynamics, and their magnitude. The same advantage was observed for onward HIV transmissions. Thus, it is known that ART decreases HIV infections [48]; hence, we can investigate the effect of ART on incidence across different age groups between men and women.

The estimates of temporal incidence were not what we expected; the difference between benchmark values and estimates values are important. If Simpact Cyan can capture estimates of age-mixing patterns in sexual partnerships, and the distribution of onward HIV transmission, the reasons for poor estimates for temporal gender and age group strata of incidence trend could be a limitation due to the summary features we used. They could not depict incidence at different time points given the dynamic nature of sexual and transmission networks. We can see that, in the Appendix B, for summary features some were not converging, as we can see in Figure A2 and Figure A3. Such instability impedes the accuracy of post-calibration. Thus, when using Simpact Cyan for real-life estimation or prediction, given available summary features to be used, an exploratory study about the accuracy level of the outputs is advisable, and parameters which need to be turned.

Between the two options of post-calibration simulation: with vectors of retained parameter values, and with medians values of retained parameter values, we produced good estimates if we ran simulations several times (1120) with the median values of the retained parameter values. However, when we checked per each parameter, the number of retained parameter values was less than 1120, which we expected with 30% of the 3752 values; this was due to the fact that sometimes the simulation may produce NAs when the model executes events which are beyond the maximum cut-off number of events (see the configuration of population.maxevents in Simpact Cyan documentation [35]). This was conducted to discard unnecessary simulations with a lot of events in order to optimise memory use. Thus, the first option of the post-calibration simulation had less simulations run compared to the second option. The more we run the simulations, more we converge to the benchmark values.

Concerning the fact that the findings of temporal trends did not meet our expectations, we assume that perhaps the summary features we used were not good enough to depict age- and gender-stratified incidence. Therefore, an exploratory simulation study is required to determine the best summary features which would better suit estimation of temporal incidence trends. Such an exploratory study could help to provide a data-driven conclusion on this matter. However, since Simpact Cyan is a multipurpose tool which can simulate many epidemiological estimates at once, even before conducting exploratory studies on summary statistics, it is crucial to undertake a sensitivity analysis [49] for specific outputs of interest before calibration.

Simpact Cyan is an agent-based model tool, which makes it computationally intense and memory expensive [50]. With simulation of viral evolution, and computation of the phylogenetic tree, the memory and time demand increases. As such, the use of combined features come at a cost in terms of time and memory [30]. Thus, the results suggest that epidemiological and sexual behavioural features were enough to estimate HIV transmission network determinants.

Therefore, if we want to estimate or predict HIV-related measurements which do not necessarily need to use the phylogenetic tree, we can use other data sources which can be captured in Simpact Cyan to avoid computational constraints. However, when we do not have other data except phylogenetic tree data, we can use features from the phylogenetic tree to calibrate the simulation framework.

We have also seen that combined features give reliable results; thus, we can use epidemiological and behavioural, and phylogenetic tree features, in methodological research, if we want to investigate accuracy of phylogenetic tree-based inferences for HIV transmission networks. The combination of epidemiological and behavioural, and phylogenetic tree features, is a bridge between two components of the simulation framework model world, namely, HIV transmission network and viral evolution.

The results of this simulation study are a note of caution for calibration simulation studies to estimate some key parameters or key determinants related to the HIV epidemic using Simpact Cyan. This study has shown that all calibration scenarios are good in terms of accuracy of estimates in post calibration analysis. However, we have seen also that not every parameter can be estimated with high accuracy given different summary features. Thus, it is advisable to first generate benchmark data and check the level of accuracy for estimates of interest using the simulation framework given the available summary features from real world data. This will help to know to what extent the results can be trusted.