Comparison of two simulators for individual based models in HIV epidemiology in a population with HSV 2 in Yaoundé (Cameroon)

Model comparisons have been widely used to guide intervention strategies to control infectious diseases. Agreement between different models is crucial for providing robust evidence for policy-makers because differences in model properties can influence their predictions. In this study, we compared models implemented by two individual-based model simulators for HIV epidemiology in a heterosexual population with Herpes simplex virus type-2 (HSV-2). For each model simulator, we constructed four models, starting from a simplified basic model and stepwise including more model complexity. For the resulting eight models, the predictions of the impact of behavioural interventions on the HIV epidemic in Yaoundé-Cameroon were compared. The results show that differences in model assumptions and model complexity can influence the size of the predicted impact of the intervention, as well as the predicted qualitative behaviour of the HIV epidemic after the intervention. These differences in predictions of an intervention were also observed for two models that agreed in their predictions of the HIV epidemic in the absence of that intervention. Without additional data, it is impossible to determine which of these two models is the most reliable. These findings highlight the importance of making more data available for the calibration and validation of epidemiological models.

In this study, the formation hazard is described by the following formula: where Pman and Pwoman are the number of partners the man and the woman already have, respectively.
The parameter represents a baseline value for relationship formation (0.1 by default). The parameters , and , represent the influence of the number of partners the man and the woman already have, respectively (0 by default). Decreasing , and , decreases the formation hazard. Therefore these parameters are used to describe the behavioural interventions in this study (see section 6 of this Supplementary Material).
To avoid that the number of relationships automatically increases with the population size, a normalization factor F is applied, which roughly divides the hazard by the size of the population.
The dissolution hazard in this study is described by a baseline value 0 for relationship dissolution (0.1 by default): ℎ = exp ( 0 ).
In this comparison study, the default settings of Simpact Cyan 1.0 for the formation and dissolution hazard were applied for the scenarios without behavioural intervention (the developers of Simpact Cyan 1.0 did not have access to the behavioural data of the city of Yaoundé described in subsection 1.2.1).
Because Simpact Cyan 1.0 keeps track of the history of each individual, we can reconstruct the sexual network from the model output. The network can be reconstructed at a given time point by including all relationships a person has at that time point for each person in the population. Furthermore, a cumulative network over a certain period of interest can be reconstructed, including all relationships between individuals in the population during that period (e.g. all relationships people had during the past 12 months). Since StepSyn 1.0 models explore behavioural differences between individuals, behavioural data on the distribution of the number of non-marital partners in the last 12 months were gathered for Yaoundé during the development of the StepSyn 1.0 modelling framework. These data were obtained from the surveys of the 4 Cities Study and provided by the Institute of Tropical Medicine (ITM) in Antwerp, Belgium (Anne Buvé, unpublished database). These data were used to estimate the parameters of a power-law distribution to generate the sexual network for the four StepSyn 1.0 models (see section "Models used in this study" in the Methods section of the main text). This work was performed before the current study.

Generation of the sexual network
Power law distributions (separately for men and women) were fitted to the data on the reported number of partners in the last 12 months taken from the 4 Cities Study. Each individual is assigned a preferred degree (number of partners within 12 months) drawn from the appropriate distribution. The individual-specific preferred degrees (PDi) are translated into probabilities of forming a new short term relationship per week (PFi), which are governed by the shortage of short links to reach the preferred degree and calculated using the formula PFi = PDi / (MD + 52), with MD being the mean duration of short term relationships in weeks, which equals 52 weeks [1]. The probability of break-up of a short link is 1/MD. With these parameters, in each week, the total male demand for new relationships is higher than the female one, partly because of the male-biased sex ratio, and partly because of female underreporting of the number of partners in the survey that serves as the basis for the distribution used. The extra male demand is uniformly distributed between the females, a method that is acceptable for modelling purposes in the absence of data on mixing patterns and life-cycle changes in sexual behaviour [2]. Although the assumption of uniformity can lead to heterogeneity in female activity being understated, assuming that the extra male demand is distributed only to the more active women would over-estimate the number of partners of females with high activity, while leaving the number of partners of females with low activity unchanged [2]. The number of weekly sex acts in married couples and short-term relationships is calculated by applying two different Poisson distributions so that their means would be consistent with the values found in the 4 Cities Study for Yaoundé [1,3].
The parameter for the proportion of male pending short links that are fulfilled by females is equal to 1 by default, which means that all male pending short links are fulfilled. Lowering this parameter decreases the probability that relationships are formed. Therefore, this parameter was used to describe the behavioural interventions (see section 6 of this Supplementary Material).

Simpact Cyan 1.0 HIV natural history
In Simpact Cyan 1.0, an infected person will go through the following HIV stages [7]: • acute HIV infection of 12 weeks [19] (parameter in Simpact Cyan 1.0: chronicstage.acutestage); • chronic HIV infectionuntil 1 year before dying of AIDS [20,21] (parameter in Simpact Cyan 1.0: agestage.start); • AIDS stage: 1 year before dying of AIDS until 6 months before dying of AIDS (default); • final AIDS stage, the person is too ill to be sexually active: 6 months before dying of AIDS until AIDS-related death (default).
In general, the formula for survival time in Simpact Cyan 1.0 is [7]: Where Vsp is a person's set point viral load. In this study x = 0 (default), k = 0 (parameter in Simpact Cyan 1.0: mortality.aids.survtime.k) and C = 2 (parameter in Simpact Cyan 1.0: mortality.aids.survtime.C). The time to AIDS-related death is implemented as the time of infection plus the survival time (here: 2 years, C = 2).  Figure S1. Interpolation spline through the data of Table 2 to obtain the HIV prevalence in 1997 for pregnant women. The smoothing parameter was determined using Generalized Cross-Validation (GCV). The figure was generated using R software version 3.6.0. (R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/). Table S1 gives an overview of the Simpact Cyan 1.0 parameters that were fitted to the calibration targets in Table 3, together with the initial range ([minimum, maximum]) used for fitting. In this formula V is the current viral load of the HIV infected person; W is a binary factor which is 1 if the uninfected person is a woman, and 0 if the uninfected person is a man; 2 is a binary factor which is 1 if the HIV infected person has HSV-2, and 0 otherwise; 2 is a binary factor which is 1 if the HIV uninfected person has HSV-2, and 0 otherwise. Table S2 gives an overview of the StepSyn 1.0 parameters that were fitted to the calibration targets in Table  3, together with the initial range ([minimum, maximum]) used for fitting.

Parameter fitting methodology
The parameters for the Simpact Cyan 1.0 models (Table S1) and the StepSyn 1.0 models (Table S2) were fitted to the calibration targets in Table 3 by applying an iterative active learning approach [8] using the procedure described in [9] and minimizing the sum of squared relative errors [10] to determine model performance. The remaining model parameters were drawn from the literature (see section 4 of this Supplementary Material).

Latin hypercube sampling
For each of the HIV transmission parameters, values were drawn from a uniform distribution, using the ranges from Table S1 and Table S2. These uniform distributions were used as initial (prior) probability distributions for the parameters in the parameter fitting procedure. We applied Latin Hypercube Sampling (LHS) [11] to select 10,000 parameter sets.

Goodness-of-fit (GOF) statistic
To find the values of the HIV transmission parameters most supported by the data in Table 3, we calculate the sum of squared relative errors [10] for each of the 10,000 parameter combinations.

Statistical evaluation of the model parameters
For the statistical evaluation of the model parameters, we focus on the subset of parameter combinations corresponding to the top 1% of the lowest values for the GOF statistic. In brief, the parameter fitting procedure consists of the following steps. First, a parameter wise comparison between the density of the initial probability distribution of the parameters and the density of the probability distribution of the subset of parameters corresponding with the top 1% solutions is conducted to determine which parameters are highly influenced by the data. Second, classification trees and generalized additive models are applied to determine which patterns of parameter vectors characterize the subspace of the top 1% solutions. Finally, we apply the Maximal Information Coefficient (MIC) [12] to determine associations between the parameters in the subspace of the top 1%.
Based on the results of the analyses above, we can narrow the solution space and repeat LHS and the steps described above several times.
A more detailed description of each step in the parameter fitting methodology, together with examples for Simpact Cyan 1.0 and StepSyn 1.0 is given below. To determine which parameters are highly influenced by the data, we compared the (prior) density of the initial uniform distribution (for the range in Tables S4 and S5) with the (posterior) density of the distribution of the top 1% parameter combinations for each parameter separately. A more peaked density for the top 1% parameter combinations indicates a parameter that is more influenced by the data.

Example 1: Simpact Cyan 1.0 basic model (Si_Ba)
If we apply univariate explorative analysis to the 10,000 LHS-generated combinations of parameters for the Si_Ba model, using the ranges in Table S1, we obtain the smoothed density plots shown in Figure S2. We observe that the baseline parameter (hiv_a) is the most influenced by the data. The remaining parameters are equally influenced by the data.  Figure S2 is a smoothed version of the histogram presented in Figure S3.

Example 2: StepSyn 1.0 basic model (St_Ba)
If we apply univariate explorative analysis to the 10,000 LHS-generated combinations of parameters for the St_Ba model, using the ranges in Table S2, we obtain the smoothed density plots shown in Figure S5. We observe that the baseline parameter (baseline) is the most influenced by the data. The remaining parameters are equally influenced by the data.

Activity region finder -top 1% parameter combinations
To determine which patterns of parameter vectors characterize the subspace of top 1% solutions, we applied Activity Region Finder (ARF) [13], a recursive partitioning classification tree algorithm. We consider each of the 10,000 parameter vectors generated by LHS as input, and the binary variable that equals 1 for a parameter combination in the top 1% and 0 otherwise as output.

Example 1: Simpact Cyan 1.0 basic model (Si_Ba)
If we applied ARF to the 10,000 LHS-generated combinations of parameters for the Si_Ba model using the ranges in Table S1, we observed that the baseline parameter (hiv_a) is responsible for the first split of the tree, followed by splits using the influence of HIV-uninfected person (exposed) being infected with HSV-2 (hiv_e2) (left and middle branch of the tree) and the influence of gender (hiv_f1) (right branch of the tree). In total 7 regions were classified as significant high activity regions (= regions associated with a low sum of squared relative errors), see Table S3. Table S3. Simpact Cyan 1.0 basic model (Si_Ba). Parameter regions classified as high activity regions by ARF. Parameters: baseline: HIV baseline probability per sex act; fm_ratio: HIV transmission probabilityinfluence of gender; hsv2_index: HIV transmission probability -influence of the HIV-infected person (index) being infected with HSV-2 (simplified HSV-2 options); hsv2_exposed: HIV transmission probability -influence of the HIV-uninfected person (exposed) being infected with HSV-2 (simplified HSV-2 options). L = left branch; M = middle branch; R = right branch. If we applied ARF to the 10,000 LHS-generated combinations of parameters for the St_Ba model using the ranges in Table S2, we observed that the baseline parameter (baseline) is responsible for the first two splits of the tree, followed by splits using the influence of HIV-uninfected person (exposed) being infected with HSV-2 (hsv2_exposed) (middle left and middle right branch of the tree) and the influence of gender (fm_ratio) (right middle branch of the tree). In total 7 regions were classified as significant high activity regions (= regions associated with a low sum of squared relative errors), see Table S4. Table S4. StepSyn 1.0 basic model (St_Ba). Parameter regions classified as high activity regions by ARF. Parameters: baseline: HIV baseline probability per sex act; fm_ratio: HIV transmission probabilityinfluence of gender; hsv2_index: HIV transmission probability -influence of the HIV-infected person (index) being infected with HSV-2 (simplified HSV-2 options); hsv2_exposed: HIV transmission probability -influence of the HIV-uninfected person (exposed) being infected with HSV-2 (simplified HSV-2 options). L = left branch; M = middle branch; R = right branch.

Generalized additive models -top 1% parameter combinations
As a second method to determine which patterns of parameter vectors characterize the subspace of top 1% solutions, we applied generalized additive models (GAM) [14] considering the same input and output variables as used with ARF.
For selecting the tuning parameter of the GAM, we consider both the Akaike Information Criterion (AIC) [15] and the Bayesian Information Criterion (BIC) [16].

Example 1: Simpact Cyan 1.0 basic model (Si_Ba)
If we applied GAM to the 10,000 LHS-generated combinations of parameters for the Si_Ba model using the ranges in Table S1, the models with tuning parameter 1 and 10 had the lowest AIC and BIC respectively, see Figure S5.    If we applied GAM to the 10,000 LHS-generated combinations of parameters for the St_Ba model using the ranges in Table S2, the models with tuning parameter 1 and 10 had the lowest AIC and BIC respectively, see Figure S8.

Maximal Information Coefficient -top 1% parameter combinations
To determine associations between the parameters within the subset of the top 1% parameter combinations, we applied the Maximal Information Coefficient (MIC) [12].  Table S5 shows the results for the top 1% of 10,000 LHS-generated combinations of parameters for the Si_Ba model using the ranges in Table S1. We observed the highest value of the MIC for the association between the baseline parameter (hiv_a) and the parameter for the influence of the HIV-uninfected person being infected with HSV-2 (hiv_e2). The second highest value of the MIC was observed for the association between the baseline parameter (hiv_a) and the parameter for the influence of the HIV-infected person being infected with HSV-2 (hiv_e1). Figure S11 shows that both associations are negative. All other associations had a MIC smaller than 0.35.

Active learning
From the results of the statistical evaluation of the parameters described in 3.3.3, we determine new ranges for the parameters and repeat the steps described in 3.3.1, 3.3.2 and 3.3.3. We repeat this several times until the solution cannot be improved anymore.
This approach is called iterative active learning [8]. In most cases, 3-4 iterations are sufficient to determine the final solution for the parameters.
In our study, the new ranges for the parameters were chosen so that they include -the parameter combination with the lowest RSSE -the significant high activity regions determined by ARF -results from GAM (only in case they define smaller intervals for some parameters than ARF)             [17] * Transmission probability from the literature converted to a transmission hazard using the procedure described in Supplementary Material, subsection 4.11.    [17] non-AIDS mortalityshape parameter for Weibull distribution mortality.normal.weibull.shape 1.8 estimated from population pyramid Yaoundé [18] non-AIDS mortalityscale parameter for Weibull distribution mortality.normal.weibull.scale 52 estimated from population pyramid Yaoundé [18] non-AIDS mortalityinfluence of genderfor a woman, half this value is added to the scale parameter of the Weibull distribution. For a man, the same amount is subtracted. mortality.normal.weibull.genderdiff 2 estimated from population pyramid Yaoundé [18] * Transmission probability from the literature converted to a transmission hazard using the procedure described in Supplementary Material, subsection 4.11.    hsv2.tr.prob.primary.m2f 0.08 estimated in this study, so that the seroprevalence of HSV-2 first increases and afterwards stabilizes at approximately 50% for females [17]   hsv2.tr.prob.primary.m2f 0.18 estimated in this study, so that the seroprevalence of HSV-2 first increases and afterwards stabilizes at approximately 50% for females [17] birth rate natality.rate 0.036 birth rate Cameroon: 3.6%* (non-AIDS) mortality rate mortality.rate 0.01 mortality rate Cameroon: 1%* immigration rate immigration.rate 0.142 popul growth in 1997: 6.8% [4]; assume emigr ~= 10%; gives immigr = 14.2% to obtain growth = 6.8% emigration rate emigration.rate 0.1 popul growth in 1997: 6.8% [4]; assume emigr ~= 10%; gives immigr = 14.2% to obtain growth = 6.8% *https://www.indexmundi.com/cameroon/demographics_profile.html  Table  3) HSV2 seed in 1980 hsv2seed.fraction 0.25 [25] * Transmission probability from the literature converted to a transmission hazard using the procedure described in Supplementary Material, subsection 4.11.  [20,22] HSV-2 transmissioninfluence of gender hsv2.tr.prob.f.m.ratio 0.5 [23,24] HIV seed in 1980 hiv.seed.strain1.initial.general 0.005 for fitting: constant seed 100 simulations: seed drawn from a binomial distribution with probability = 0.005 leads to an HIV prevalence of 0.6% for males and 1.2% for females in 1989 (values in Table 3) HSV2 seed in 1980 hsv2.seed.initial.general 0.25 [25]

Converting transmission probability parameters to hazard parameters
While in the majority of the literature, and also in StepSyn 1.0, transmission parameters are described as probabilities, the parameters of Simpact Cyan 1.0 are described in terms of hazards. Transmission probability parameters were converted to hazard parameters using the following formula [26]: where F(t) is the cumulative distribution function and λ(x) is the hazard function.
An example of conversion of a transmission probability parameter to a hazard parameter is presented below.  Table 3.  Figure S13 shows the median HSV-2 for females and the range of 100 simulations for the period 1980-2005 for the 8 models.

Behavioural interventions
For the Simpact Cyan 1.0 models, the parameters describing the weight for the number of relationships a person already has ( , and , in formula (1) (parameter name in Simpact Cyan 1.0: formation.hazard.agegap.numrel_man, formation.hazard.agegap.numrel_woman) was changed. For the StepSyn 1.0 models, we changed the parameter for the proportion of male pending short links (concurrent unstable relationships) that are fulfilled by females to reach the individual preferred degree (parameter name in StepSyn 1.0: pending.short.links.fulfilled)(see Supplementary Material, section 1.2.2. for more detail).
Because sexual networks for Simpact Cyan 1.0 and StepSyn 1.0 are generated in different ways, the initial sexual network before applying interventions was different between Simpact Cyan 1.0 and StepSyn 1.0 models. Figure S13 shows the distribution of the number of partners at the start of the HIV epidemic for both Simpact Cyan 1.0 and StepSyn 1.0 in case no intervention was implemented. In Simpact Cyan 1.0 the parameter for the weight for the number of relationships a person already has is equal to 0. In StepSyn 1.0 the parameter for the proportion of male pending short links that are fulfilled by females is equal to 1.   Changing the weight for the number of relationships a person already has from 0 to -0.05 in Simpact Cyan 1.0, and changing the proportion of male pending short links that are fulfilled by females from 1 to 0.7 in StepSyn 1.0 reduces the mean number of partners by 6% and the 95th percentile with 1, while keeping the median and the 75th percentile unchanged.
For all models, behavioural interventions also reduced the cumulative HSV-2 incidence (see Table S22).        Table S24. Results for the predicted HIV prevalence in 2004 for females and males in case of no behavioural intervention, a behavioural intervention implemented in 1990, and a lower promiscuity from 1980 onwards. Interventions were implemented as follows. For Simpact Cyan 1.0, the weight for the number of relationships a person already has (α_(numrel,man) = α_(numrel,woman) in formula (1)) was changed from 0 to -0.05. For StepSyn 1.0, the proportion of male pending shorts links fulfilled by females (pending.short.links.fulfilled) was changed from 1 to 0.7. Median HIV prevalence (in %) and range ([minimum,maximum]) of 100 simulations. Models abbreviations as in Table 4 of the main text.

Lifetime number of partners
Estimates from the literature for the lifetime number of partners in Yaoundé were compared with outcomes from the model simulations. Ferry et al. 1  The reported information in Ferry et al. 1 was based on a questionnaire from UNAIDS 28 , which determined the lifetime number of partners as the number of partners up to the date participants filled in the questionnaire.
If we assume that individuals are equally spread across ages 15-49 years, and that people become sexually active at the age of 15 years, then the mean period of sexual activity for the population in Ferry et al. and the 2004 DHS is 17 years.
The models without inflow and outflow assume a population that is sexually active over de whole simulation period 1980-2005. So the total number of partners during the whole simulation is based on a period of sexual activity of 25 years (25 years of sexual activity is simulated for all individuals). We therefore multiply the total number of partners for the models with no inflow and outflow with 17/25.
For the models with inflow and outflow, the period of sexual activity that is modeled varies between 0 and 25 year (depending on when the person enters or leaves the population). This means that the mean period of sexual activity is 12.5 years if we assume that individuals are equally spread across ages. We correct for this by multiplying the total number of partners in these models by 17/12.5. Table S26. Lifetime number of partners reported in the literature and estimated from the models (after correcting for period of sexual activity). Models abbreviations as in Figure 1 of the main text.