Understanding how Victoria, Australia gained control of its second COVID-19 wave

During 2020, Victoria was the Australian state hardest hit by COVID-19, but was successful in controlling its second wave through aggressive policy interventions. We calibrated a detailed compartmental model of Victoria’s second wave to multiple geographically-structured epidemic time-series indicators. We achieved a good fit overall and for individual health services through a combination of time-varying processes, including case detection, population mobility, school closures, physical distancing and face covering usage. Estimates of the risk of death in those aged ≥75 and of hospitalisation were higher than international estimates, reflecting concentration of cases in high-risk settings. We estimated significant effects for each of the calibrated time-varying processes, with estimates for the individual-level effect of physical distancing of 37.4% (95%CrI 7.2−56.4%) and of face coverings of 45.9% (95%CrI 32.9−55.6%). That the multi-faceted interventions led to the dramatic reversal in the epidemic trajectory is supported by our results, with face coverings likely particularly important.


Platform for infectious disease dynamics simulation
We developed a deterministic compartmental model of COVID-19 transmission using the AuTuMN platform, publicly available at https://github.com/monash-emu/AuTuMN/. Our repository allows for the rapid and robust creation and stratification of models of infectious disease epidemiology and includes pluggable modules to simulate heterogeneous population mixing, demographic processes, multiple circulating pathogen strains, repeated stratification and other dynamics relevant to infectious disease transmission. The platform was created to simulate TB dynamics, being an infectious disease whose epidemiology differs markedly by setting, such that considerable flexibility is desirable [1]. We have progressively developed the structures of our platform over recent years, and further adapted it to be sufficiently flexible to permit simulation of other infectious diseases, such as COVID-19. A similar model has been applied to several countries of the Asia-Pacific, with the application to the Philippines (without structure to represent geospatial stratification or contact tracing) previously described. [2] Base COVID-19 model Using the base framework of an SEIR model (susceptible, exposed, infectious, removed), we split the exposed and infectious compartments into two sequential compartments each (SEEIIR). The two sequential exposed compartments represent the non-infectious and infectious phases of the incubation period, with the latter representing the "presymptomatic" phase such that infectiousness occurs during three of the six sequential phases. For this reason, "active" is a more accurate term for the two sequential "I" compartments and is preferred henceforward. The two infectious compartments represent early and late phases of active disease, during which symptoms occur if the disease episode is symptomatic, and allow explicit representation of notification, case isolation, hospitalisation and admission to the intensive care unit (ICU). The "active" compartment also includes some persons who remain asymptomatic throughout their disease episode, such that these compartments do not map directly to either persons who are infectious or those who are symptomatic ( Figure 1).
The latently infected and infectious presymptomatic periods together comprise the incubation period, with the incubation period and the proportion of this period for which patients are infectious defined by input parameters described below. In general, two sequential compartments can be S E E I I R used to form a gamma-distributed profile of transition to infectiousness following exposure if the progression rates for these two compartments are equal, although in implementing this model the relative sojourn times in the two sequential compartments usually differed. Nevertheless, the profiles implemented are broadly consistent with the empirically observed log-normal distribution of individual incubation periods [3].
The transition from early active to late active represents the point at which patients are detected (for those persons for whom detection does eventually occur) and isolation then occurs from this point forward (i.e. applies during the late active phase only). This transition point is also used to represent the point of admission to hospital or transition from hospital ward to intensive care for patients for whom this occurs (see Section 1).

Age stratification
All compartments of this base compartmental structure were stratified by age into five-year bands from 0-4 years of age through to 70-74 years of age, with the final age group being those aged 75 years and older. Heterogeneous baseline contact patterns by age were incorporated using agespecific contact rates estimated from survey data from the POLYMOD study, which reported rates of age-specific contacts in various locations between persons of different age groupings (as described below 1). These are then modified by non-pharmaceutical interventions. Our modelled age groups were chosen to match these mixing matrices. The automatic demographic features of AuTuMN that can be used to simulate births, ageing and deaths were not implemented, because the issues considered pertain to the short-to medium-term and the immediate implementation of control strategies, for which population demographics are less relevant.

Contact matrices construction
For each location L (home, school, work, other locations) the age-specific contact matrix C L = (c L i,j ) ∈ R 16×16 + is defined such that c L i,j is the average number of contacts that a typical individual aged i has with individuals aged j. As there is no contact survey available for Australia that is complete across all age groups, the matrices C L were obtained by extrapolating contact matrices from the United Kingdom, being a country included in the POLYMOD study in 2005 [4]. The original matrices from the United Kingdom are denoted Q L = (q L i,j ) ∈ R 16×16 + , where q L i,j is defined using the same convention as for c L i,j . The matrices Q L were extracted using the R package "socialmixr" (v 0.1.8) and then adjusted to account for age distribution differences between Victoria and the United Kingdom.
Let π j denote the proportion of people aged j in Victoria, and ρ j the proportion of people aged j in the United Kingdom. The contact matrices C L were obtained from: In sensitivity analysis, we replaced the United Kingdom with Belgium as the proxy country from which survey data were used for constructing Victoria's matrices. Belgium was chosen as the POLYMOD country with the closest similarity in age distribution of Australia and Belgium in 2005.

Clinical stratification
The age-stratified late exposed/incubation and both the early and late active disease compartments were further stratified into five "clinical" categories: 1) asymptomatic, 2) symptomatic ambulatory, never detected, 3) symptomatic ambulatory, ever detected, 4) ever hospitalised, never critical and 5) ever critically unwell ( Figure 2). The proportion of new infectious persons entering stratum 1 (asymptomatic) is age-dependent. The proportion of symptomatic patients (strata 2 to 5) ever detected (strata 3 to 5) is set through a parameter that represents the time-varying proportion of all symptomatic patients who are ever detected (the case detection rate). Of those ever symptomatic (strata 2 to 5), a time-constant but age-specific proportion is considered to be hospitalised (entering strata 4 or 5). Of those hospitalised (entering strata 4 or 5), a fixed proportion was considered to be critically unwell (entering stratum 5, Figure 3).   to ICU admission.

Infectiousness
Asymptomatic persons are assumed to be less infectious per unit time active than symptomatic persons not undergoing case isolation (typically by around 50%, although this is varied in calibration/uncertainty analysis). Infectiousness is also decreased for persons who have been detected to reflect case isolation, and for those admitted to hospital or ICU to reflect infection control procedures (by 80% for both groups). Presymptomatic individuals are assumed to have equivalent infectiousness to those with early active COVID-19.

Application of COVID-19-related death
Age-specific infection fatality rates (IFRs) were applied and distributed across strata 4 and 5, with no deaths typically applied to the first three strata. A ceiling of 50% is set on the proportion of those admitted to ICU (entering stratum 5) who die. If the infection fatality rate is greater than this ceiling, the proportion of critically unwell persons dying was set to 50%, with the remainder of the infection fatality rate then applied to the hospitalised proportion. Otherwise, if the infection fatality rate is less than half of the absolute proportion of persons critically unwell, the infection fatality rate is applied entirely through stratum 5 (such that the proportion of critically unwell persons dying in that age group becomes <50% and the proportion of stratum 4 dying is set to zero). In the event that the infection fatality rate for an age group is greater than the total proportion hospitalised (which is unusual, but could occur for the oldest age group under certain parameter configurations), the remaining deaths are assigned to the asymptomatic stratum. This approach was adopted for computational ease and is valid because the duration active for persons entering this stratum is the same as for the other non-hospitalised strata, such that the dynamics are identical to assigning the deaths to any of the first three strata. We used the age-specific IFRs previously estimated from age-specific death data from 45 countries and results from national-level seroprevalence surveys [5]. Uncertainty in the IFR estimates used in the model is incorporated as described in Section 10.

Determining the proportion of cases detected
We calculate a time-varying case detection rate, being the proportion of all symptomatic cases (clinical strata 2 to 5) that are detected (clinical strata 3 to 5). This proportion is informed by the number of tests performed using the following formula: time is the time in days from the 31 st December 2019 and tests(time) is the number of tests per capita done on that date. To determine the value of the shape parameter, we solve this equation based on the assumption that a certain daily testing rate tests(t) is associated with a certain CDR(t). Solving for shape yields: That is, if it is assumed that a certain daily per capita testing rate is associated with a certain proportion of symptomatic cases detected, we can determine shape. As this relationship is not well understood and unlikely to be consistent across all settings, we vary the CDR that is associated with a certain per capita testing rate during uncertainty/calibration. Given that the CDR value can be varied widely, the purpose of this is to incorporate changes in the case detection rate that reflect the empirical historical profile of changes in testing capacity over time.
The proportion of persons entering clinical stratum 3 is calculated once the CDR is known, along with the proportion of all incident cases hospitalised (strata 4 and 5).

Isolation of detected cases
As described in the Section 1 above, as infected persons progress from the early to the late stage of active COVID-19, infectiousness is reduced for those in the detected strata (3 to 5) to reflect case isolation.

Model adaptation
We simulate quarantining of persons identified as first degree contacts of COVID-19 patients explicitly through stratification of the compartments representing active COVID-19. That is, the compartments representing both phases of the incubation period and both phases of active COVID-19 are duplicated, with model strata referred to as "traced" and "untraced". In model initialisation, all infectious seed is assigned to the untraced stratum. All newly infected persons commence their incubation period in the untraced stratum of the early incubation period. As for isolated and hospitalised patients, those undergoing quarantine have their infectiousness reduced by 80%.

Contact tracing process
Identification of infected persons through contact tracing is assumed to apply to those in their early The proportion of contacts of identified cases that is traced, q(t), is considered to decrease as the severity of the COVID-19 epidemic increases, because we expect contact tracing to decline in efficiency as more cases are identified. That is, we assume that contact tracing is universal as COVID-19 prevalence approaches zero and declines exponentially with increasing prevalence. The relationship between the proportion of contacts of identified patients who are quarantined and prevalence is given as: Rather than estimate τ directly, we estimate the more intuitive quantity of the proportion of contacts of identified patients who would be quarantined at a particular prevalence. Solving for the previous equation for τ , we obtain: at a specific prevalence that accords with a particular value of q. Fixing prev 0 at 10 −3 , we can vary q 0 in calibration as the proportion of contacts of identified cases detected at a prevalence of one active case per thousand population.
Finally, q(t) × u(t) gives the proportion of all infected persons who are traced. This proportion of persons entering their early latent period transition to the equivalent compartment in the traced stratum before proceeding to the late latent period. Figure 4 presents the time-varying processes and parameters relating to contact tracing implementation in the model.

Testing data
Statewide daily testing data by date of test were provided by DHHS and applied to all health service clusters to provide a broad profile of the variation in testing capacity over time, including the lower testing numbers in early June compared to at the peak of the epidemic. Data sparseness precluded us from implementing separate functions for each individual health service cluster. For this application to Victoria, the case detection proportion corresponding to a per capita rate of testing of one test per thousand population per day was varied as a calibration parameter in creating the time-varying case detection proportion function. Note that testing rates were typically considerably higher than one per thousand per day during the period modelled, such that the actual modelled case detection proportion is considerably higher than the case detection calibration parameter for most of the simulation period.

Implementation of non-pharmaceutical interventions
A major part of the rationale for the development of this model was to capture the past impact of non-pharmaceutical interventions (NPIs) and produce future scenarios projections with the implementation or release of such interventions.

Isolation and quarantine
For persons who are identified with symptomatic disease and enter clinical stratum 3, self-isolation is assumed to occur and their infectiousness is modified as described above. The proportion of ambulatory symptomatic persons effectively identified through the public health response by any means is determined by the case detection rate as described above.

Community quarantine or "lockdown" measures
For all NPIs relating to reduction of human mobility or "lockdown" (i.e. all NPIs other than isolation and quarantine), these interventions are implemented through dynamic adjustments to the age-assortative mixing matrix.
The age-specific mixing matrices we describe above (Section 1) have the major advantage of allowing for disaggregation of total contact rates by location, i.e. home, work, school and other locations. This disaggregation allows for the simulation of various NPIs in the local context by dynamically varying the contribution of each location to reflect the historical implementation of the interventions.
The corresponding mixing matrix (denoted C 0 ) is presented using the standard convention that a row represents the average number of age-specific contacts per day for a contact recipient of a given age-group. In other words, the element C 0i,j is the average number of contacts per day that an individual of age-group j makes with individuals of age-group i .
This matrix results from the summation of the four location-specific contact matrices: where C H , C S , C W and C L are the age-specific contact matrices associated with households, schools, workplaces and other locations, respectively.
In our model, the contributions of the matrices C S , C W and C L vary with time such that the input contact matrix can be written: The modifying functions are each squared to capture the effect of the mobility changes on both the infector and the infectee in any given interaction that could potentially result in transmission.
The modifying functions incorporate both macro-distancing and microdistancing effects, depending on the location.

School closures/re-openings
Reduced attendance at schools is represented through the function s(t), which represents the proportion of all school students currently attending on-site teaching. If all education is fully remote, s(t) = 0 and C S does not contribute to the overall mixing matrix C(t). s(t) is calculated through a series of estimates of the proportion of students attending onsite education, to which a smoothed step function is fitted. Note that the dramatic changes in this contribution to the mixing matrix with school closures/re-openings is a more marked change than is seen with the simulation of policy changes in workplaces and other locations (which are determined by empiric data and so do not vary so abruptly or reach a value of zero).

Workplace closures
Workplace closures are represented by quadratically reducing the contribution of workplace contacts to the total mixing matrix over time. This is achieved through the scaling term w(t) 2 which modifies the contribution of C W to the overall mixing matrix C(t). The profile of the function w(t) is set by fitting a polynomial spline function to Google mobility data for workplace attendance (Table 2).

Community-wide movement restriction
Community-wide movement restriction (or "lockdown") measures are represented by proportionally reducing the contribution of the other locations contacts to the total mixing matrix over time.
This is achieved through the scaling term l(t) 2 which modifies the contribution of C L to the overall mixing matrix C(t). The profile of the function l(t) is set by fitting a polynomial spline function to an average of Google mobility data for various locations, as indicated in Table 2.

Household contacts
The contribution of household contacts to the overall mixing matrix C(t) is fixed over time. Although Google provides mobility estimates for residential contacts, the nature of these data are different from those for each of the other Google mobility types in that they represent the time spent in that location rather than the number of attendances. The daily frequency with which people attend their residence is likely to be close to one and we considered that household members likely have a daily opportunity for infection with each other household member. Therefore, we did not implement a function to scale the contribution of household contacts to the mixing matrix with time.

Microdistancing
Interventions other than those that prevent people coming into contact with one another are thought to be important to COVID-19 transmission and epidemiology, such as maintaining interpersonal physical distance and the wearing of face coverings. We therefore implemented a "microdistancing" function to represent reductions in the rate of effective contact that is not attributable to persons visiting specific locations and so is not captured through Google mobility data. This microdistancing function reduces the values of all elements of the mixing matrices by a certain proportion. These time-varying functions multiplicatively scale the location-specific contact rate modifiers s(t), w(t) and l(t).

School closures
The effect of Victorian school closures is captured through the timeline presented in Table 3.  Table 3: Timeline used to implement Victorian school closure policies. The function is applied to both metropolitan and regional services.

Macrodistancing in workplaces and other locations
The functions applied here are determined by the Google mobility data according to Table 2, as described above, but are applied separately for each service. Because Google mobility data pertains to local government areas (LGAs), whereas health service clusters may receive patients from across the state, it was necessary to map mobility data to services. Health service clusters' overall mobility values in each location were calculated using a weighted average of LGA mobility values according to the historical pattern of the origin of patients presenting to services within each service.
As a hypothetical example, if 50% of patients historically presenting to Barwon South West health services come from the City of Geelong, the mobility data for the City of Geelong will contribute 50% of the Google mobility estimate of Barwon South West.
Historical patterns of patient presentations by health service cluster were provided by the Victorian Department of Health and Human Services (DHHS).

Microdistancing approach
In this application to Victoria, the microdistancing function m(t) is comprised of two components: physical distancing and face coverings. Both physical distancing and face coverings microdistancing are applied to the three non-household locations, such that the microdistancing function for non-household locations is given by: The two interventions are assumed to be independent and so are multiplicative. As for the macrodistancing functions, the two functions of time are squared to represent their effects on both the infector and the infectee in any potentially infectious interaction.

Physical distancing
The physical distancing function d(t) is a transposed and translated hyperbolic tan function. The parameters of this function were estimated by using maximum a posteriori inference, with priors that penalised large shape parameters (to avoid extremely rapid transitions). The proportions of respondents answering "always" to YouGov surveys of Victorian residents asking "Thinking about the last 7 days, about how many people from your household have you come into physical contact with (within 2 meters / 6 feet)?" were used as input data. Resulting parameters were: shape, 0.262764; lower asymptote, 0.2803973; upper asymptote; 0.4421819; and inflection point, 15 th July. The resulting function is presented in Figure 5.

Face coverings
Two separate face coverings microdistancing functions are employed, one for metropolitan and one for regional health service clusters. These functions were fitted using the same methods as for physical distancing, using YouGov data on Victorian residents' survey responses to the question "Thinking about the last 7 days, have you worn a face mask outside your home (e.g. when on public transport, going to a supermarket, going to a main road)?". Estimated parameters were: shape, 0.5261693; lower asymptote, 0.130469; upper asymptote, 0.9143849; and inflection point, 26 th July (consistent with the policy change in metropolitan Melbourne). This was applied directly to metropolitan services and translated ten days later for regional services, where face coverings were mandated from the 2 nd August. The resulting function is presented in Figure 6.

Between service mixing
The preceding section describes the creation of heterogeneous mixing matrices by age for each of the nine health service clusters individually. These mixing matrices are then combined to create a single time-varying heterogeneous mixing matrix by service and age resulting in a 144 by 144 (9 × 16 = 144) square mixing matrix. The force of infection for an index service is calculated from the mixing matrices of the age-assortative matrix for each of the services modelled. The spatial mixing matrix is based on the adjacency of health service clusters as indicated in Table 4.

Model initialisation
The model was commenced from around one to two weeks earlier than the actual beginning of Victoria's second wave (as determined by genomic analysis), in order that the distribution of infectious persons distributes naturally across compartments as the model approaches the actual beginning of Victoria's second wave in early June. The actual start date selected was the 14 th May. The infectious seed needed at this time was then calibrated to ensure dynamics were realistic at the beginning of the second wave. The infectious seed is distributed evenly across metropolitan services, consistent with the epidemic's emergence from Metropolitan Melbourne.  [6,7,8,9]. A systematic review [3] found that data are best fitted by a log-normal distribution (mean 5.8 days, CI 5.0 to 6.7, median 5.1 days). Our systematic review [10] found that estimates of the mean incubation period have varied from 3.6 to 7.4 days.

Non-age-stratified parameters
continuation of parameters

Age-specific parameters
Age-structured parameters are presented in Table 6.

Incidence
Incidence is calculated as any transitions into the early active compartment ("I").

Hospital occupancy
This is calculated as the sum of three quantities: 1. All persons in the late active compartment in clinical stratum 4, representing those admitted to hospital but never critically unwell.
2. All persons in the late active compartment in clinical stratum 5, representing those currently admitted to ICU.
3. A proportion of the early active compartment in clinical stratum 5, representing those who will be admitted to ICU at a time in the future. This proportion is calculated as the quotient of 1) the difference between the pre-ICU period and the pre-hospital period for clinical stratum 4, divided by 2) the total pre-ICU period. That is, a proportion of the pre-ICU period is considered to represent patients in hospital who have not yet been admitted to ICU.

ICU occupancy
This is calculated as all persons in the late active compartment in clinical stratum 5.

Seropositive proportion
This is calculated as the proportion of the population in the recovered ("R") compartment. Although very similar numerically to the attack rate, persons who died of COVID-19 are not included in the denominator.

COVID-19-related mortality
This is calculated as all transitions representing death, exiting the model. This is implemented as depletion of the late active compartment.

Notifications
Local case notifications are calculated as transitions from the early to the late active compartment for clinical strata 3 to 5.

Calibration
We calibrated the model using the adaptive Metropolis algorithm described by Haario et al. [18].
A standard Metropolis algorithm with fixed proposal distribution parameters was used for the first 500 iterations to initiate the covariance matrix before the adaptive algorithm commenced. Seven chains were then run to ensure 10,000 iterations post-burn-in were achieved.

Rationale for service-specific targets
For all services (both metropolitan and regional), we included the time series of daily notifications for that service as a calibration target, using a normal distribution for the likelihood function. A normal distribution is preferred because the mapping process for the notifications for each service results in these quantities not being integer-valued.
In addition, we include time series for the following quantities at the state level.

Variation of age-specific proportion hospitalised using "adjuster" parameters
Our parameters included the age-specific proportions of cases hospitalised, which were varied up and down together during calibration. These proportion parameters were modified through an "adjuster" parameter that is not strictly a multiplier, but is rather implemented in such a way as to scale the base parameter value while ensuring that the adjusted parameter remains a proportion (with range zero to one). This adjuster parameter can be considered as a multiplicative factor that is applied to the odds ratio that is equivalent to the baseline proportion to be adjusted. Specifically, the adjusted proportion is equal to: proportion × adjuster proportion × (adjuster − 1) + 1

Variation of infection fatality rate
The infection fatality rate can be defined as the risk of death given an episode of infection, including asymptomatic and undetected episodes. This is considered a more stable quantity than the case fatality rate. However, it is still likely to vary considerably between settings and so is adjusted during the calibration process. Because the epidemic in Victoria was characterised by high rates of transmission and disease in aged care, we fix the infection fatality rate for all age groups up to 74 years to the estimates derived from the literature, but vary the infection fatality rate for those aged 75 and above. This is intended to capture the increased risk of death for those in residential aged care during the second wave, the large majority of whom would be included in this age bracket.

Likelihood function
Likelihood functions are derived from comparing model outputs to target data at each time point nominated for calibration.
The composite likelihood function is given formally as: where t indexes the date, g indexes the service, n t refers to daily new notifications, d t to daily deaths, h t to daily new hospitalisations and i t to daily new ICU admissions. The contributions of each state-wide component to the composite likelihood are measured with Poisson distributions (e.g. n t (θ) = P oiss(ν t (θ)), where ν t (θ) is the number of notifications simulated by the model at date t under parameter set θ), and normal distributions are used for each n t,g (because these targets are not integer-valued). σ is the ratio of the peak of each service-specific notification to the corresponding standard deviation of each of the normal distributions used in calculating their contribution to the likelihood. This was included as a calibration parameter to improve calibration efficiency.

Ordinary differential equations
For the clearest description of the model, we refer the reader to our code repository, because our object-oriented approach to software development is intended to be highly transparent and readable. For those who prefer dynamical systems such as this presented in the form of ordinary differential equations, we present the following. dS a,g dt = −λ a,g (t)σ a S a,g dE a,g,q=1 dt = λ a,g (t)σ a S a,g − αE a,g,q=1 − χ(t)E a,g,q=1 dE a,g,q=2 dt = −αE a,g,q=2 + χ(t)E a,g,q=1 dP a,c,g,q dt = p a,c (t)αE a,g,q − νP a,c,g,q dI a,c,g,q dt = νP a,c,g,q − γ c I a,c,g,q dL a,c,g,q dt = γ c I a,c,g,q − δ a,c L a,c,g,q − µ a,c L a,c,g,q re g (t) + gr g (t) + pa g (t) + tr g (t) 4 Symbol Explanation

Sensitivity analysis using alternative base mixing matrices
We ran a sensitivity analysis using age-specific mixing matrices derived from survey data for Belgium, rather than the United Kingdom, which were weighted to the Victorian age structure using the same methods as described above. The results of this analysis were similar to that of the baseline analysis and are presented in Figure 21, Figure 22 and Figure 23. These results show no significant epidemiological differences from those derived from the base case analysis.

Sensitivity analysis with Google residential mobility used to scale home location contribution to mixing matrices
We ran an alternative analysis, in which the home location contribution to the dynamic mixing matrices scaled with Google residential mobility data. This replaced the baseline assumption that  Table 4 the rates of home location contacts remained fixed throughout the simulations. The results of this analysis were similar to that of the baseline analysis and are presented in Figure 24, Figure 25 and With this approach, the posterior distributions for the incubation period and the duration active shortened well into the lower tail of their respective normal prior distributions. Meanwhile, the effect of face coverings was clustered towards the upper bound of its prior distribution. Exploring the model through manual calibration, the reason for this was determined to be that it was not possible to achieve a sufficiently steep decline in case rates without shortening the serial interval and emphasising the effect of face coverings. This remained the case even if the effect of face coverings was allowed to reach extremely high (and implausible) levels (e.g. 100%). This is because the importance of home contacts increased towards the peak of the epidemic, as Google residential mobility increased. We do not consider this realistic because these contacts would likely have

Additional analysis of calibration convergence
In the main analysis described above, each of the seven chains of the adaptive Metropolis algorithm was initialised by adding random noise around parameter values obtained from a previous calibration. This approach was used to assist the algorithm to reach the parameter space's high posterior areas within a reasonable number of iterations. The high number of dimensions (N = 18) made other approaches such as Latin-Hypercube-based initialisation impractical, as chains could take millions of iterations before adequately sampling plausible parameter sets.
In a separate experiment, we ran another adaptive Metropolis algorithm where only ten parameters were varied, but this time initialisation was performed with Latin Hypercube Sampling (LHS) across ten independent chains. The ten parameters were selected based on their expected influence on the main findings and whether they were among the primary estimands of interest  Figure 28). We found that the chains converged to similar values when considering the LHSbased calibration compared to the original analysis. The posterior ranges obtained with the LHSbased analysis were sometimes slightly narrower compared to the original calibration. This is attributable to the fact that more parameters were fixed throughout calibration in the LHS-based calibration than in the original calibration. Consequently, more constraints were imposed on the varied parameters that were previously found to be correlated with parameters that have been