Study design

Validation cohorts

For assay validation, 338 negative control samples were obtained from arbovirus-unexposed adult French (n = 249)⁴⁶ and Peruvian donors (n = 89). Suspected alphavirus positive control samples from adult Senegalese donors (n = 38) with serological evidence of prior CHIKV infection were provided by the Centre for Arboviruses and Haemorrhagic Fever Viruses (CRORA) in Institut Pasteur de Dakar. Suspected flavivirus positive control samples (n = 222) were obtained from adult Peruvian donors (n = 102) with either suspected or confirmed DENV exposure (i.e. based on clinical history, serological testing, or PCR verification) and from adult Senegalese donors (n = 120, provided by CRORA) with serological evidence of suspected prior DENV (n = 60), WNV (n = 30), or YFV exposure (n = 30).

Cross-sectional surveys

Samples were obtained from eight cross-sectional surveys conducted at sites in Senegal, Peru, French Guiana, and New Caledonia. In Peru, 880 samples were obtained from an all-age cross-sectional population survey of arboviruses and malaria conducted in 2018 in three peri-urban communities in Iquitos district, Department of Loreto, in the Peruvian Amazon⁴⁷. In Senegal, 1219 samples were obtained from four population-level cross-sectional malaria prevalence surveys conducted in 2016 and 2018 in villages Dielmo and Ndiop, which have been monitored longitudinally since the 1990’s with regard to the transmission of malaria⁴⁸. In French Guiana, 787 samples were obtained from two multicentric cross-sectional surveys conducted in the Maroni area along the French Guiana-Suriname border in 2015 (n = 409) and in 2019 (n = 378) to study the prevalence of malaria in adult goldminers⁴⁹. In New Caledonia, 600 samples were obtained from adults 18–64 years of age participating in the 2021–2022 Adult Health Barometer, a cross-sectional descriptive survey⁵⁰.

Ethics declaration

No new samples were collected as part of this study. This study involved the secondary use of samples previously collected in other studies. For all studies, the ethical approval and the informed consent covered this additional analysis.

The surveys in French Guiana were approved by the Commissie voor Mensgebonden Wetenschappelijk Onderzoek (CMWO) in Suriname (Opinion Number VG10-14 & DVG-738) and by the Institut national de la santé et de la recherche médicale (INSERM) Ethics Evaluation Committee in France (No. 14–187 of 09.12.2014).

The Peruvian samples were collected from participants in the study Circles of Research on Arboviruses and Malaria (SIDISI 101645/2017), following approval by the Institutional Ethics Committee of Universidad Peruana Cayetano Heredia (UPCH).

The Senegalese samples were collected as part of the ongoing Dielmo and Ndiop project, which was examined and approved by the Senegalese National Health Research Ethics Committee (CNERS Sénégal). Approval to measure antibodies to arboviruses in these samples was granted by CNERS Sénégal (N 00000007 MSAS/CNERS/Sec). The study in New Caledonia was approved by the “Comité de Protection des Personnes Ouest IV” (21.05.16.92458-avis_50_21_2 and avis_50-21_2MS1) and by the Consultative Ethics Committee of New Caledonia (Avis 2021-03 002). The study was recorded on Clinicaltrials.gov (ID: NCT05218304). The negative French samples were obtained from the COVID-Oise study, which was registered with ClinicalTrials.gov (NCT04644159) and received ethical approval from the Comité de Protection des Personnes Nord Ouest IV. The negative Peruvian samples were collected from healthy individuals from Lima as part of the Molecular and Sero-epidemiology of relapsing Plasmodium vivax in Peru (SIDISI 100873/2017) study, which was approved by the Institutional Ethics Committee of UPCH.

Bead-based immunoassay procedures

A bead-based immunoassay was developed and optimised for simultaneous detection of IgG antibody responses to 28 recombinant antigens, including both structural and non-structural proteins, from commercially available sources. For CHIKV, antigens included a virus-like particle (VLP) (containing envelope, pre-membrane and membrane proteins) and the envelope protein 2 (E2). For MAYV, antigens included the E2 and the non-structural protein 1 (NSP1), whereas for ONNV, only a VLP was included.

DENV antigens included the envelope protein domain III (EDIII), a VLP, and the non-structural protein 1 (NS1) from each of the four serotypes (DENV1–4). ZIKV antigens included the EDIII and a VLP, as well as two distinct versions of NS1, representing the Suriname (SU) and Ugandan (U) viral strains, to account for geographical differences in strain-specific exposure. For YFV and JEV, antigens included NS1 and the envelope protein (E), whereas for USUV, only NS1 was included. For WNV, antigens included NS1 and EDIII. All antigens were commercially available and were produced by The Native Antigen Company, except for WNV EDIII, and MAYV NSP1, which were produced by Sino Biological and Alpha Diagnostic, respectively. Each recombinant antigen was covalently coupled to a spectrally unique set of carboxylated paramagnetic nanobeads (MagPlex; Luminex Corp) according to the manufacturer’s instructions. Optimal protein concentrations were determined experimentally to generate a log-linear standard curve with a positive control plasma pool⁵¹. All antigens, their catalogue numbers, and the protein coupling concentrations are listed in Supplementary Table 2.

The multiplex assay was performed according to a previously described protocol⁵¹. Briefly, 50 μl of each plasma sample, diluted 1:100 in phosphate-buffered saline with 1% (w/v) bovine serum albumin and 0.05% (v/v) Tween 20 (PBT buffer), was added to each well of a 96-well plate and incubated with 50 μl of bead mix (containing 250 beads per analyte per well) for 30 min. Following three washes with PBT buffer, 100 μl of the secondary antibody (R-Phycoerytrin (PE) conjugated F(ab’)₂ fragment Donkey anti-Human IgG, Jackson ImmunoResearch: JIR709116098), diluted 1:120 in PBT-buffer, was added to each well and incubated for 15 min. After three final washes with PBT buffer, beads were re-suspended in 100 μl of PBT buffer and plates were read using a Luminex MagPix instrument (Luminex Corp.). The same 11 technical controls, including pools of sera from both arbovirus-unexposed donors and exposed donors, were assayed on each plate and used to control for plate-to-plate-related assay variation.

Competitive immunoassays

Samples were pre-incubated with soluble antigens (MAYV E2, CHIKV E2, or BSA) at increasing concentrations to competitively block or deplete antibodies prior to detection on bead-bound antigen. Each sample was assayed under 30 individual conditions (3 competitor antigens × 5 concentrations × 2 bead-bound antigens). Fifty microliters of the sample at a 1:100 dilution in PBS-T was pre-incubated with each soluble competitor antigen (MAYV E2, CHIKV E2, or BSA) at 0, 0.01, 0.1, 1, and 10 µg/ml with shaking for 1 h at room temperature. After pre-incubation, each 50 µl aliquot was mixed with 50 µl of antigen-coupled beads (either MAYV E2 or CHIKV E2). Antibody responses towards MAYV and CHIKV were quantified independently in separate wells (mono-plex) but otherwise as described above (see Bead-based immunoassay procedures). MAYV E2 antibody binding was measured in the presence of soluble MAYV E2, CHIKV E2, or BSA, and CHIKV E2 binding was measured separately under the same competition conditions. The obtained MFI of each sample and experimental condition was expressed as a proportion of the maximum sample MFI, defined as the MFI measured without soluble competitor antigen.

Data management, normalisation and analysis

Data management and analysis were performed in R (R version 4.3.3). Bayesian finite mixture models were executed from R on the Stan platform (Stan Development Team, 2024, Stan Modelling Language version 2.35, https://mc-stan.org) using the cmdstanr package (cmdstanr: R Interface to ‘CmdStan’, version 0.7.1) as the R interface to Stan. Mathematical models of antibody cross-reactivity and virus transmission were implemented in Python (v. 3.13.1) using the PyMC module (v. 5.20.0).

Normalisation of immunoassay data

To account for plate-to-plate assay variation, a robust linear model was fitted to the data for each antigen from the 11 technical controls:

where \(\log \left({{Ab}}_{{ij}}\right)\) is the observed signal in logarithmic scale for technical control i on plate j; \({\beta }_{1}({\text{control}}_{i})\) represents the effect of technical control i sample accounting for differences in antibody quantity in different controls; \({\beta }_{2}({\text{plate}}_{j})\) accounts for overall differences between plates, such as differences due to protein coupling batch, or experimental and/or reading conditions; and \({\epsilon }_{{ij}}=\text{Normal}\left(0,{\sigma }^{2}\right)\) is the random error with an assumed normal distribution centred around zero with homogeneous variance across all plates. Once the best fit parameters are identified, the MFI of all samples and controls were corrected by:

where \(\log ({\text{Ab}}_{{ij}}^{{\prime} })\) is the log normalised MFI of sample i on plate j on the logarithmic scale.

Bayesian mixture models

A two-part GMM with shared parameters was fitted jointly to the log-transformed antibody data from the known negative (unexposed) control samples and all cross-sectional survey samples in a Bayesian framework. The negative control samples were assumed to be described by the lowest Gaussian component only, while the cross-sectional survey samples were assumed to be described by all the components such that:

where \(k\) represents the number of mixture components, \({\text{Ab}}_{\text{neg}}\) and \({\text{Ab}}_{\text{cs}}\) represent the antibody response in the negative control samples and the cross-sectional survey samples, respectively, and \({\theta }_{1}\),…, \({\theta }_{k}\) are the mixing coefficients (or weights) for each of the \(k\) components.

It is not uncommon that more than two subpopulations may be present in the data. In this case, a simple two-component model will not capture the complexity of the data, and therefore both two- (\(k=2\)) and three-component (\(k=3\)) mixtures were evaluated for all antibody responses. The lower component distribution and upper component distribution were assumed to represent the seronegative and seropositive populations, respectively. However, in the case of a three-component mixture, the class membership (i.e. seronegative or seropositive) of the middle component is not explicitly assigned and was determined by comparing the reactivity of the model-estimated middle component with the observed reactivity of the suspected positive control samples (to which the GMM was not fitted).

All Bayesian GMM models were fitted with four parallel Markov Chain Monte Carlo (MCMC) chains using the No-U-Turn-Sampler (NUTS) algorithm. Uniform priors were assumed for all parameters, but parameters \({\mu }_{1},\ldots,{\mu }_{k}\) were constrained so that \({\mu }_{1} . Model convergence and appropriate chain-mixing were evaluated by examining the model diagnostic plots and by assessing the R-hat and the bulk effective sample size statistics of convergence and efficiency, respectively. The best fit model for each antibody response (i.e. \(k=2\) vs \(k=3\)), to be retained for downstream analysis, was identified by visually inspecting the model fits to the data, and more formally using efficient approximate leave-one-out cross-validation (LOO) to identify the model with the lowest LOOIC⁵².

Finite mixture models based on the more flexible Weibull probability distribution were also evaluated (see Supplementary Information: “Methods” section).

ROC analysis

For each antibody response, ROC analysis was applied to assess the trade-off between sensitivity and specificity at all possible seropositivity thresholds.

An “empirical” ROC curve was calculated from the observed data from the negative and positive control samples. For the empirical ROC curve, the assay specificity is estimated as the proportion of negative control samples that are below a specified threshold, and the assay sensitivity is estimated as the proportion of positive control samples that are greater than that threshold.

In a situation where known seropositive and seronegative samples are not available, the ROC curve cannot be calculated directly from the data. It is, however, possible to generate an “analytical” ROC curve based on the parameter estimates from a finite mixture model and the CDF of the component normal distributions. For the analytical ROC curve, the specificity is estimated as the proportion of the CDF(s) of the mixture component(s) representing the seronegative population that is below a specified threshold, and assay sensitivity is estimated as the proportion of the CDF(s) of the mixture component(s) representing the seropositive population that is greater than that threshold. The performance of each antibody assay in discriminating seronegative from seropositive individuals was assessed using the AUROC, and an AUROC > 0.9 was considered acceptable.

Graphical LASSO partial correlation network analysis

Correlation network analysis was performed on log-transformed antibody responses. For each antigen, responses were first residualised using linear models with site as a predictor. Conditional-dependence networks were estimated from the residualised data using EBIC graphical LASSO based on Spearman’s rank correlation. Network layouts were computed from the resulting partial correlation matrix using the Fruchterman–Reingold force-directed algorithm. Network stability was assessed by non-parametric bootstrapping (1000 iterations), and stable edges were defined by a bootstrap inclusion probability ≥0.5 with confidence intervals excluding zero. For visual clarity, non-stable edges and edges with partial correlations 

PCA

PCA was applied to log-transformed, mean-centred antibody data for all flavivirus NS1 antigens. PC scores were extracted and visualised for the first three PCs to assess clustering of samples in reduced-dimensional space. Antibody response loadings were extracted for the same components to examine the relative contribution from each antigen. For visualisation purposes, PCA plots were colour-coded according to serological response phenotypes derived from the univariate GMM. Samples not exceeding the GMM-derived threshold of seropositivity for any of the NS1 antigens were designated “Unexposed-like”. For each NS1 antigen, samples with a posterior mean probability > 0.5 of belonging to the highest mixture component were classified as high responders and serological response phenotypes were defined as follows: High non-DENV/non-ZIKV-flavivirus, high for at least one NS1 antigen from WNV, USUV, JEV, or YFV and not high for DENV and ZIKV; High DENV-only, high for at least one DENV NS1 antigen and not high for any other; High ZIKV, high for at least one ZIKV NS1 antigen and not high for any other; High DENV + ZIKV only, High for at least one DENV and at least one ZIKV NS1 antigen and not high for any other; Pan-Flavivirus-high, high for at least one DENV and at least one ZIKV NS1 antigen and high for all other tested NS1 antigens (WNV, USUV, JEV, and YFV).

Mathematical model of antibody dynamics, cross-reactivity and virus transmission

We adapted a previously validated model developed by Hozé et al.¹² to jointly characterise antibody cross-reactivity and reconstruct individual infection histories for MAYV and CHIKV within a unified inference framework. Specifically, the model quantifies how infection with one virus shapes the antibody response to that virus and how this response modulates the measured reactivity to a related, non-infecting virus. Simultaneously, it estimates the probability that individuals fall into one of four prior infection states: uninfected, infected with MAYV, infected with CHIKV, or infected with both. We fit the model jointly to log₁₀-MFI data for MAYV E2 and CHIKV E2 from cross-sectional survey participants and known unexposed controls, where the latter contribute information on the background reactivity in unexposed individuals. We incorporate a hierarchical structure to account for site-specific variation in this background reactivity. For the subset of samples evaluated using the competitive immunoassay (n = 30), we assign a fixed infection status based on the assay results (see above).

Model overview

Let \(M\) and \(C\) refer to MAYV and CHIKV, respectively. For an individual \(j\), \(({t}_{j}^{\text{M}},{t}_{j}^{\text{C}})\) denote the antibody levels in log₁₀-MFI for MAYV and CHIKV, respectively. The latent infection status is denoted \({I}_{j}=({i}_{\text{M}},{i}_{\text{C}}) \, \epsilon \, {\left\{\mathrm{0,1}\right\}}^{2}\), where 1 indicates a past infection. \(\theta\) denotes the full parameter vector, and \(\varphi ({x|}\mu,\sigma )\) represents the normal probability density function with mean \(\mu\) and standard deviation σ.

The joint probability of the measured log₁₀-MFI and the latent infection status is given by:

Where the first term models the probability of the antibody response intensities (measured log₁₀-MFI) conditional on the infection history, and the second term models the probability of the infection status conditional on the virus transmission process, i.e. the FOI.

Antibody response component

We model the measured antibody response intensities (log₁₀-MFI) conditional on an individual’s infection history.

Site-specific background:

True unexposed individuals will still exhibit non-zero assay reactivity. The distribution of this non-specific background reactivity will vary between different populations. Therefore, the background log₁₀-MFI values are modelled hierarchically to share information across sites and from known unexposed controls. The Peruvian unexposed controls from Lima were assumed to share the same background reactivity as the Peruvian survey participants from Iquitos.

\({\mu }_{0,s}^{V}\) denotes the site-specific background log₁₀-MFI mean for virus \(V \, \epsilon \, \left\{M,C\right\}\) and is modelled hierarchically:

Boosts and cross-reactivity:

The conditional distributions of antibody responses are defined for each possible infection status as:

1. 1.

   Uninfected \((\mathrm{0,0})\)

   $$P({t}_{j}^{\text{M}},{t}_{j}^{\text{C}}|{I}_{j}=(0,0),\theta )=\varphi ({t}_{j}^{\text{M}}|{\mu }_{0}^{\text{M}},{\sigma }^{\text{M}})\times \varphi ({t}_{j}^{\text{C}}|{\mu }_{0}^{\text{C}},{\sigma }^{\text{C}})$$

   (6)

2. 2.

   Infected by MAYV only \((\mathrm{1,0})\):

   $$P({t}_{j}^{\text{M}},{t}_{j}^{\text{C}}|{I}_{j}=(1,0),\theta )=\varphi ({t}_{j}^{\text{M}}|{\mu }_{0}^{\text{M}}+{\mu }^{\text{M}},{\sigma }^{\text{M}})\times \varphi ({t}_{j}^{\text{C}}|{\mu }_{0}^{\text{C}}+{\gamma }^{\text{M}\to \text{C}}{t}_{j}^{\text{M}},{\sigma }^{\text{C}})$$

   (7)

3. 3.

   Infected by CHIKV only\((\mathrm{0,1})\):

   $$P({t}_{j}^{\text{M}},{t}_{j}^{\text{C}}|{I}_{j}=(0,1),\theta )=\varphi ({t}_{j}^{\text{M}}|{\mu }_{0}^{\text{M}}+{\gamma }^{\text{C}\to \text{M}}{t}_{j}^{\text{C}},{\sigma }^{\text{M}})\times \varphi ({t}_{j}^{\text{C}}|{\mu }_{0}^{\text{C}}+{\mu }^{\text{C}},{\sigma }^{\text{C}})$$

   (8)

4. 4.

   Infected by both MAYV and CHIKV \((\mathrm{1,1})\):

Here \({\mu }^{\text{M}}\), \({\mu }^{\text{C}}\) are homologous additive boosts to the log₁₀-MFI signal following infection and \({\gamma }^{\text{M}\to \text{C}}\) and \({\gamma }^{\text{C}\to \text{M}}\) are multiplicative cross-reactive effects whereby infection with one virus increases the observed log₁₀-MFI against the heterologous virus.

Virus circulation component

\({\lambda }_{s}^{V}(y)\) denote the annual FOI for virus \(V\) at site \(s\) in year \(y\). For an individual \(j\) of age \({a}_{j}\), residing in site \({s}_{j}\), the cumulative hazard due to virus \(V\), experienced up to the survey year \(S\), is given by:

FOI models:

We consider two models for the FOI time series:

1. 1.

   Constant FOI model:

   $${\lambda }_{s}^{V}(y)={\lambda }_{s}^{V}$$

   (11)

2. 2.

   Outbreak FOI model:

where:

\({T}_{s}^{V}\) is the year of peak transmission, \({\tau }_{s}^{V}\) controls the temporal width of the outbreak, and \({\alpha }_{s}^{V}\) is the total FOI associated with the outbreak.

Epidemiological situation:

Based on the available epidemiological evidence, we assumed a single CHIKV outbreak to have occurred in 2014 in Maroni (French Guiana)³⁹ and in New Caledonia⁴¹. We assumed MAYV transmission to be constant over time in Maroni¹² and CHIKV transmission to be constant over time in Dielmo and Ndiop (Senegal)³⁸. We assumed MAYV to be absent from Senegal (Dielmo and Ndiop) and New Caledonia. For Peru (Iquitos), we fitted alternative transmission models for MAYV and CHIKV, assuming either a constant FOI or a single outbreak. The most likely scenario was identified based on model fit using the WAIC, alongside inspection of MCMC chain convergence.

Infection probability

Assuming independence between the infection processes of the two viruses, the joint infection status \({I}_{j}=({i}_{M},{i}_{C})\) probabilities follow the standard serocatalytic model. For example, the probability of being infected by MAYV only is:

Likelihood

The full likelihood comprises three data components:

1. 1.

   Unknown infection status:

   For each individual \(j\epsilon E\) with unknown infection status, the likelihood contribution is summed over all possible infection statuses \({I}_{k}\epsilon \left\{(\mathrm{0,0}),(\mathrm{1,0}),(\mathrm{0,1}),(\mathrm{1,1})\right\}\):

   $${L}_{exp }=\mathop{\prod }\limits_{j \, \epsilon \, E}\mathop{\sum }\limits_{{I}_{k}}P({t}_{j}^{\text{M}},{t}_{j}^{\text{C}}|{I}_{j}={I}_{k},\theta )\times P({I}_{j}={I}_{k}|\theta )$$

   (14)

2. 2.

   Assay-confirmed infection status:

   For individuals \(j \, \epsilon \, A\) with infection status confirmed by the competitive assay \({I}_{j}^{{assay}}\), the contribution is:

   $${L}_{\text{assay}}=\mathop{\prod }\limits_{j \, \epsilon \, A}P({I}_{j}={I}_{j}^{\text{assay}}|\theta )$$

   (15)

3. 3.

   Unexposed controls:

For individuals \(j \, \epsilon \, U\) known to be uninfected, the likelihood is:

The full log-likelihood is then

Alternative model incorporating heterologous cross-protection

To assess whether prior infection with one virus modifies the risk of acquiring the other, we extended the standard serocatalytic framework by introducing multiplicative cross-protection parameters. \({\lambda }_{s}^{\text{M}}(y)\) and \({\lambda }_{s}^{\text{C}}(y)\) denote the annual FOI for MAYV and CHIKV, respectively, at site \(s\) and year \(y\). If an individual has previously been infected with CHIKV, their subsequent hazard of MAYV infection is reduced to \({\lambda }_{s}^{\text{M}}\times (1-{\delta }^{\text{C}\to \text{M}})\), where \({0\le \delta }^{\text{C}\to \text{M}}\le 1\) quantifies the protection conferred by CHIKV antibodies. Similarly, a prior MAYV infection reduces the subsequent CHIKV infection hazard to \({\lambda }_{s}^{\text{C}}\times (1-{\delta }^{\text{M}\to \text{C}})\).

We assumed that co-infection within the same year is not possible, and both FOIs remain constant over one-year intervals. Under these assumptions, the joint infection probabilities are updated in discrete one-year time steps using the modified hazards:

where \({P}_{({i}_{\text{M}},{i}_{\text{C}})}(t)\) denotes the probability of being in infection state \(({i}_{\text{M}},{i}_{\text{C}})\) at time t.

Model comparison using WAIC^53,54 showed no improvement in fit from including heterologous cross-protection (\({\varDelta }_{{WAIC}}\le -1.6\)), so we retained the simpler model for inference.

Inference

Parameters were estimated in a Bayesian framework using MCMC using the NUTS algorithm implemented in the PyMC module⁵⁵. Convergence was assessed through visual inspection of trace plots.

Priors

\({T}_{s}^{\text{C}}\) Fixed to 2014 for Peru, New Caledonia, and French Guiana (single outbreak model).

\({\tau }_{s}^{\text{C}}\) Fixed to \(\frac{1}{2}\) so that 95% of the outbreak FOI lies within ~2 years.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.