Modeling lesion transition dynamics to clinically characterize mpox patients in the Democratic Republic of the Congo

Description of cohort and study design

We used comprehensive clinical data collected from MPXV-infected patients at the remote L’Hopital General de Reference de Kole (Kole Hospital), in the rainforest of the Congo River basin of DRC, from March 2007 until August 2011. The data were collected as part of a study conducted jointly by the Institute National de Recherche Biomedical (INRB) and the US Army Medical Research Institute of Infectious Diseases (USAMRIID) [27].

Of 244 patients diagnosed with MPXV infection, 228 tested positive by MPXV-specific PCR testing. In addition, 13 patients were excluded because of incomplete observational data. In our analysis, we studied lesion transition dynamics in 149 participants and virus dynamics in 151 participants, excluding 66 and 64 individuals, respectively, owing to incomplete longitudinal data on lesion counts and viral load (see Extended Data Fig 1). Among the 228 patients who tested positive, 3 were fetal cases. The dataset comprised lesion counts at different stages on various areas of the patients’ bodies, as well as the MPXV genomic DNA load (i.e., viral load) in peripheral blood, pharyngeal swabs, and scabs from the same individuals. Additionally, we integrated clinical laboratory data from hematology analysis, urinalysis, and clinical chemistry within the same patient cohort; Table 1 shows the summary statistics of these measured items for all individuals (total) and by the level of disease severity (G1, G2, and G3), which is characterized in the next section. Furthermore, the comprehensive clinical dataset for these patients included annotations for their sex, age, clinical symptoms, clinical signs, and exposure possibilities (details are outlined in [27]).

• Download figure
• Open in new tab

Extended Data Figure 1 Flowchart of the cohort, along with the number of participants and inclusion criteria for our analysis.

Quantifying and stratifying lesion transition dynamics of MPXV-infected patients

Generally, mpox lesions can be classified into seven distinct stages: macules, papules, vesicles, pustules, umbilication, scabbing, and desquamation [27, 28, 32]. As depicted in the distribution of lesion counts at various stages in Extended Data Fig 2, the average lesion count at the macule, papule, vesicle, umbilicated, and pustule stages reached zero by approximately 20 days after onset of the first lesion, whereas the count at the scab and desquamation stages persisted even beyond 40 days post-onset (Fig 1A). Nevertheless, there was significant individual-level variation in the number of lesions and in the temporal dynamics of lesion transition (Extended Data Fig 2).

• Download figure
• Open in new tab

Extended Data Figure 2 Time course distribution of the average lesion count of each stage:

Time course plots of individual lesion counts for each stage (i.e., macule, papule, vesicle, pustule, umbilicated, scab, and desquamation) across all patients.

• Download figure
• Open in new tab

Figure 1 Stratification and characterization of lesion transition dynamics:

(A) Time course plot of the average lesion count for each stage (i.e., macule, papule, vesicle, pustule, umbilicated, scab, and desquamation) across all patients. (B) UMAP of stratified lesion transition dynamics based on the estimated parameters. Data points represent individual participants and are colored according to group (i.e., G1, G2, and G3: green, blue, and red). (C) The time-course patterns highlighted by Partial Least-Squares Discriminant Analysis (PLS-DA), which discriminates each group from the others. All lesion count values for each time were reconstructed to highlight the difference between lesion stages. The mean (solid line) and mean ± 1 standard deviation (shaded area) calculated based on these reconstructed values are shown, colored according to stage (i.e., L₁, L₂, L_3, and L₄: orange, magenta, light purple, and cyan). (D) Total number of lesions at each stage for each stratified group, calculated based on the reconstructed lesion transition dynamics. (E) Number of lesions (left) and period from lesion onset to disappearance (right), calculated based on reconstructed lesion dynamics. For (D) and (E), statistical significance was calculated using the pairwise Mann-Whitney U test. Also, p-values were corrected by Bonferroni’s method (NS.: p-value > 0.05, ∗: p-value ≤ 0.05, ∗∗: p-value ≤ 0.01, and ∗∗∗: p-value ≤ 0.001, respectively). The statistical significance of the difference is indicated on the top of the plots. (F) Kaplan-Meier survival curves of lesion disappearance for each group based on reconstructed lesion dynamics. Solid lines and shaded areas indicate the survival probabilities and 95% confidence intervals, respectively. The green, blue, red, and black colors correspond to G1, G2, G3, and the total group, respectively. The table below the figure shows the number and percentage of survivors in each group at each time point. (G) Confusion matrices on the relationship between the stratified group and CDC severity score (top panel, individuals with <25, 25-99, 100-499, and >=500 lesions) as well as between stratified group and age category (bottom panel, individuals <5, 5-11, and >=12 years old).

To quantify the variation among MPXV-infected patients (e.g., Extended Data Fig 2), we used a mathematical model that describes the dynamics of lesion transition (see Methods in detail). For model fitting, we merged the macule, papule, and vesicle stages into one stage, owing to the limited number of lesions observed at each stage and the lack of clinical relevance in distinguishing between these three stages. Additionally, considering that the infectious period persists until all scabs (or crusts) have fallen off [28], we excluded the desquamation stage from our analysis. Consequently, we modelled lesion counts using the compartments L₁, L₂, L₃, and L₄, in which L₁ represents lesions in the macule, papule, and vesicle stages, and L₂, L_3, and L₄ represent lesions in the pustule, umbilicated, and scab stages, respectively (Extended Data Fig 3). We then characterized the transition dynamics among these stages at the individual level. We reconstructed the lesion transition dynamics for each participant in our cohort over a period of approximately 40 days following the onset of lesions (see Methods for details). The individual-specific model fits are shown in Fig S1 and the estimated parameters are summarized in Table S1.

• Download figure
• Open in new tab

Extended Data Figure 3 Schematic diagram of compartment model for lesion transition:

The lesion counts in the macule, papule, and vesicle (i.e., these three are in one compartment, as described in the main text); pustule; umbilicated; and scab stages are defined as L₁, L₂, L₃, and L₄, respectively. It is assumed that L₁ increases after t₀ days post lesion onset according to a per-lesion rate ge^−ct as the infection progresses (with L₁(t) = 1 for t ≤ t₀), where g corresponds to the maximum growth rate of the number of lesions and c determines the rate at which this growth rate decreases over time. In contrast, L₄ decreases according to a rate d, and the transition rates between the stages are denoted as a_n(1 − e^−τ_nt) for n = 1,2,3 where a_n is the maximum transition rate and 1 − e^−τ_nt sets the time dependence of the delay between the stage transitions.

We performed an unsupervised clustering analysis using a dissimilarity-based random forest clustering approach [37, 38] to stratify the time-course patterns of lesion transition based on estimated individual-level parameters (see Methods in detail). The number of clusters optimized by the algorithm was three (i.e., three groups of patients; G1, G2, and G3). Fig 1B represents a two-dimensional Uniform Manifold Approximation and Projection (UMAP) embedding of these three groups, which clearly shows that G1 is distinct from the other groups. Using a different color for each group, we also plotted the reconstructed individual transition dynamics in Extended Data Fig 4. The corresponding projected time course of each group in Fig 1C shows that the lesion transitions between stages progressed rapidly among individuals in G1 (N=42), whereas they were slower in the G3 participants (N=31). On the other hand, individuals in G2 (N=76) showed intermediate patterns, even though the dynamics of the scab disappearance stage (the cyan curve) were comparable to that in G3. These patterns suggest substantial variation in the infectious period of MPXV-infected patients, consistent with our recent findings of substantial individual variations in the duration of viral shedding in different specimens of MPXV clade IIb B.1 [39]. Our statistical tests showed that individuals in G1 and G3 had significantly smaller and larger numbers of lesions for all stages, respectively, compared with individuals in the other groups (Fig 1D). Additionally, we observed statistically significant differences in total lesion counts and lesion duration (defined as time between onset and lesion disappearance) between patient groups, except between G2 and G3 in lesion duration (Fig 1E). These findings demonstrate that G2 and G3 are distinguished by the total number of lesions as well as by the lesion counts in each stage, rather than by lesion duration. Hereafter we classify G1 as “mild progressors” and G3 as “severe progressors” in our analysis.

• Download figure
• Open in new tab

Extended Data Figure 4 The time-course patterns of reconstructed lesion dynamics for each stage:

The reconstructed individual lesion dynamics for each stage, shown for each stratified group (i.e., G1, G2, and G3).

To describe the difference in timing of lesion disappearance, we calculated the probability of detectable lesions over time by using the model with estimated parameters for each group (Fig 1F). We found that individuals in G1 had significantly shorter periods of lesion symptoms: there was a difference in lesion disappearance between G1 vs G2, and between G1 vs G3 (p<0.0001 by log-rank test for both). In the total group (i.e., a group of all analyzed cases), the probability of detectable lesions was 55% (95% CI: 48%-64%) at 31 days after lesion onset, consistent with a reported median duration of experiencing apparent lesions among mpox cases in the 2022 outbreak [40]. All lesions in G1 patients had disappeared at 31 days after lesion onset, whereas the corresponding probabilities in G2 and G3 were 51% (95% CI: 41%-64%) and 71% (95% CI: 57%-89%), respectively. In all stratified groups, the probability of detectable lesions was less than 10% by 44 days after lesion onset.

Additionally, we compared our stratified groups with two major categorizations of MPXV-infected patients to assess reliability from a clinical standpoint: "CDC [Centers for Disease Control and Prevention] lesion severity categories” [25] based on the total number of lesions on admission day (i.e., individuals with <25, 25-99, 100-499, and >=500 lesions are categorized as 1, 2, 3, and 4, respectively) and "age categories” [27] (individuals <5, 5-11, and >=12 years old are categorized as 1, 2, and 3, respectively) (Fig 1G). Both categories were based on clinical investigations. We observed that CDC severity categories generally aligned with our stratified groups (e.g., patients in G1 typically had fewer lesions at admission, and those in G3 had more lesions). In our logistic regression analysis, using the number of lesions at symptom onset as the predictor and the stratified group as the outcome, we obtained high ROC-AUC values: 94% for G1 and 82% for G3. Maximizing the Youden Index from the ROC curves, we calculated cutoff values of 107 and 281 lesion counts with 62% and 61% specificity and 87% and 72% sensitivity for predicting mild and severe progressors, respectively. Notably, the cutoff value for G1 (i.e., 107 lesion counts) closely matched the boundary between CDC lesion severity categories 2 and 3 (100 lesion counts) [28]. In terms of age categories, children under 5 years old generally exhibited a longer lesion duration and a larger number of lesions, which aligns with previous reports indicating severe lesion symptoms in young children [17, 27, 41].

Viral load at lesion onset may define severity on lesion symptoms

While the CDC lesion severity categories are useful, it is practically (and ethically) challenging in clinical settings to always count all lesions appearing on the entire body of a patient [42]. Recent studies have suggested that the diagnosis of small lesions is difficult even for clinicians and quite subjective [43]. To provide more objective support for the mpox diagnostics, we explored the relationship between the lesion transition dynamics and individual-level viral load, as our recent study suggested that viral load in lesions is positively associated with the infectious period of a patient (i.e., the duration of having culturable MPXV) for clade IIb MPXV [33, 34, 44].

In our cohort, although we measured the viral load in peripheral blood, pharyngeal swabs, and scabs [27] (see Methods), our focus was on the viral load in peripheral blood because the other viral loads are not routinely measured (typically once per patient in most cases). All patients (N=151; Extended Data Fig 1) showed a decaying viral load profile (Fig 2A), since the viral load is likely to have already peaked before the first measurement (i.e., before patients present to the hospital) [27]. Therefore, we used a decay model to reconstruct the best-fit viral load profile for each case, taking into account inter-patient heterogeneity (see Methods). In Fig 2B, we depict the reconstructed individual dynamics aligned with the stratified groups (i.e., G1, G2, and G3). We compared four key features derived from the viral dynamics among the groups in Fig 2C: viral load at lesion onset and at disappearance (when scab count falls below 1), total viral load (defined as the AUC of viral load above 10⁻² genomes/mL), and the viral clearance rate. Interestingly, we found that the onset viral loads and total viral load in G1 were significantly lower than those in G3, whereas the viral loads at lesion disappearance in G1 were significantly higher than those in G2. Conversely, there was no apparent difference in the clearance rate among the groups. On an individual patient level, we observed positive correlations for the two main quantified lesion dynamics features (lesion duration and total lesion counts in Fig 1E) with the onset viral load in Fig 2D (0.39 and 0.49 Pearson’s correlation coefficient, respectively, with a p-value of 0.001 in both cases). These findings suggest that the onset viral load may serve as a potential biomarker for predicting lesion transition patterns and severity of illness over the course of infection (discussed below). Additional correlation analyses between lesion dynamics and viral load dynamics are illustrated in Extended Data Fig 5.

• Download figure
• Open in new tab

Extended Data Figure 5 Relationship between lesion dynamics and virus dynamics:

Correlations between the features of individual-level lesion transition dynamics (i.e., the lesion duration and total lesion counts) and those of individual-level decay viral load dynamics (i.e., the viral clearance rate, the total viral load and the viral load at lesion disappearance).

• Download figure
• Open in new tab

Figure 2 Relationship between lesion transition dynamics and peripheral blood viral load:

(A) The measured individual decay dynamics of mpox viral load (i.e., genomic DNA) in peripheral blood samples obtained from the DRC cohort. (B) The individual-level mpox decay viral load, reconstructed using the mathematical model, plotted aligned with the stratified groups (i.e., G1, G2, and G3: green, blue, and red). (C) Comparison of viral load at onset, viral clearance rate, total viral load in peripheral blood, and viral load at lesion disappearance among the stratified groups. (D) Correlations between the features of individual-level lesion transition dynamics (i.e., the lesion duration and total lesion counts) and those of individual-level decay viral load dynamics (i.e., the viral load at lesion onset and disappearance). (E)(F) The ROC curves for a generalized linear model (GLM) logistic regression with the viral load at lesion onset for predicting stratified group G1 or not and G3 or not, respectively. The corresponding ROC-AUCs are displayed in the bottom right of each panel. FPR and TPR mean false positive rate and true positive rate, respectively. (G) Swarm plot of the viral load at onset among the stratified groups. The dashed line indicates the threshold value (4.6 log₁₀ genomes/mL). (H) Comparison of lesion duration between two groups classified based on the threshold value (Low: below 4.6 log₁₀ genomes/mL, High: above 4.6 log₁₀ genomes/mL).

Finally, we explored whether an individual’s risk group (G1, G2, or G3) could be predicted by the onset viral load in peripheral blood. We used a logistic regression analysis with the onset viral load as the predictor variable and the stratified groups as the outcome variable. Although we could not attain a high ROC-AUC for predicting G2 (41%), which exhibits intermediate patterns, we did achieve relatively higher ROC-AUCs of 69% for G1 (defined as mild progressors) and 67% for G3 (defined as severe progressors) as shown in Fig 2EF. Similarly, based on the Youden Index, we calculated cutoff values for the onset viral load of 4.557 and 4.619 (log₁₀ genomes/mL) with 61.8% and 61.1% specificity and 87.1% and 72.0% sensitivity for predicting mild and severe progressors, respectively. Interestingly, these two cutoff values were very close, and therefore we henceforth use a common threshold of an onset viral load of 4.6 (log₁₀ genomes/mL) instead. Fig 2G shows that the threshold generally distinguishes between G1 and non-G1 (i.e., mild or non-mild progressors) as well as between G3 and non-G3 (i.e., severe or non-severe progressors), but not the intermediate pattern (i.e., G2). Notably, the threshold viral load significantly differentiated the duration of lesions (p<0.0001 by Mann-Whitney U-test): 28.0 and 35.6 days on average are required until lesion disappearance for patients with an onset viral load below and above 4.6 (log₁₀ genomes /mL), respectively (Fig 2H).

Characterizing mild and severe progressors by clinical laboratory data

In this section, we further explored potential biomarkers for use in predicting mild and severe progressors by incorporating additional clinical laboratory data. Using the clinical laboratory data obtained from our cohort [27] (see Methods), we annotated the 149 patients with hematology analysis, urinalysis, and clinical chemistry on the day of lesion onset (Table 1 and Extended Data Fig 1). We examined the correlation between the stratified group based on estimated lesion transition patterns and the 41 variables (summarized in Table 1). Each factor was compared between the three groups by Kruskal-Wallis and Fisher’s exact test where the p-values were corrected by the false discovery rate (FDR). However, we found that only the white blood cell count differed significantly among the stratified groups (i.e., adjusted p-value of less than 0.05).

Next, we investigated whether the inclusion of the clinical laboratory data alongside the onset viral load improved the performance of the model predictions in Fig 2EF. We here used a supervised machine learning algorithm (Light Gradient Boosting Machine [LightGBM]) to predict the group from the laboratory data of 120 patients in our cohort (see details in Methods). Despite the inclusion of additional data, we achieved similar ROC-AUCs, specifically for G1 (70%) and G3 (64%), as depicted in (Fig 3AB), along with a lower ROC-AUC for G2 (53%). Interestingly, the combination of viral load data and clinical laboratory information did not significantly improve predictability for either mild or severe progressors.

• Download figure
• Open in new tab

Figure 3 Clinical characterization of stratified groups:

(A) and (B) depict the ROC curves of LightGBM classifiers trained to predict either G1 (i.e., mild progressor) or not, and G3 (i.e., severe progressor) or not, respectively. The corresponding ROC-AUCs are calculated and displayed in the bottom right of each panel. (C) and (D) show bar plots of how often a feature was selected among seeds as an indicator of selection stability. Selected features are shown with a frequency greater than 0.25. OnsetVL: viral load at lesion onset, MCHC: mean corpuscular hemoglobin concentration, GLU: blood glucose level, UA: uric acid, PMN: neutrophils, MCH: mean corpuscular hemoglobin, MCV: mean corpuscular volume, AST: aspartate aminotransferase, LYM: lymphocytes, MONO: monocytes, WBC: white blood cell, EOSIN: eosinophils, GGT: gamma glutamyl transferase.

Finally, concentrating on G1 (mild progressors) and G3 (severe progressors), we examined features (other than the onset viral load) that predict which of these two non-intermediate groups’ individuals belong to. In Fig 3AB, we selected features for our model to enhance classification accuracy and efficiency by eliminating redundancy. Specifically, the feature selection was based on the importance of SHapley Additive exPlanations (SHAP) values [45], which were added to the LightGBM model in order of significance (see Methods). The features with the highest ROC-AUC in the training dataset were then chosen. Following feature selection, the frequencies of the features selected over 50 iterations [46] for G1 and G3 are shown in Fig 3CD. Additionally, the distribution of corresponding features between G1 as opposed to non-G1, and G3 as opposed to non-G3, is depicted in Extended Data Fig 6, respectively. The frequency lists of selected features for both G1 and G3 included the onset viral load as the 1st and 8th feature, respectively. This again supports the importance of onset viral load in predicting the severity of lesion symptoms. Regarding other biomarkers, we found that low viral load and low AST (which indicates absence of liver damage at onset) were weakly associated with mild progressors. In contrast, for the severe progressors, the most relevant biomarkers were high WBC and AST, which implies elevated inflammatory response and liver dysfunction at onset. These clinical findings suggest that inflammatory response and virus-induced organ damage at onset may be linked to disease progression.

• Download figure
• Open in new tab

Extended Data Figure 6 Distribution of features selected by machine learning:

(A) and (B) show distributions of selected features by LightGBM with G1 or not, and G3 or not, respectively. Statistical significance is calculated using the pairwise Mann-Whitney U test. (NS.: p-value > 0.05, ∗: p-value ≤ 0.05, ∗∗: p-value ≤ 0.01, and ∗∗∗: p-value ≤ 0.001, respectively). The statistical significance of the difference is indicated on the top of the plots.

Ethics statement

The prospective observational study of the clinical natural history of human MPXV infections at the remote Kole Hospital in the rainforest of the Congo River basin of the DRC, which was conducted from March 2007 until August 2011, was reviewed and approved by the Human Use Committee of the United States Army Medical Research Institute of Infectious Diseases (FY05-13) and the Headquarters, United States Army Medical Research and Development Command Institutional Review Board (IRB), Frederick, MD, USA, as well as the Ethics Committee at the University of Kinshasa School of Public Health (KSPH), Kinshasa, DRC (see [27] for details). The present study was approved by the ethics committee of Nagoya University (approval number: Hc 22-07).

Study data

The count of lesions on nine distinct skin regions and the oropharynx of each patient’s body (a total of 10 areas) was documented and reported in [27]. In this study, we used the total number of lesions for each stage (macules, papules, vesicles, pustules, umbilication, scabbing, and desquamation) by aggregating the lesion count from all areas for each day of the study. Additionally, we incorporated viral load data, measuring the amount of MPXV genomic DNA in peripheral blood, pharyngeal swabs, and scabs. Clinical laboratory data, obtained through hematology analysis, urinalysis, and clinical chemistry from the same group of patients, was also considered. Detailed information on data preparation can be found in [27].

Modeling lesion transition dynamics

In our recent report [27], we categorized lesions into seven different stages: macule, papule, vesicle, pustule, umbilicated, scab, and desquamation. Because of the limited number of macule, papule, and vesicle lesions for the purpose of parameter estimation and the lack of clinical relevance in distinguishing between these three stages, we merged these into one stage. The notation L₁(t) therefore represents the total number of macule, papule, and vesicle lesions at time t since lesion onset (i.e., t = 0 represents the day at which a patient reported its first lesion appearance). Additionally, considering that the infectious period persists until all scabs or crusts have fallen off [28], we considered lesion counts in the umbilicated, pustule, and scab stages as L₂(t), L₃(t), and L₄(t), respectively, while excluding lesion data from the desquamation stage. Subsequently, we developed the following mathematical model to describe the dynamics of lesion transitions in MPXV-infected patients:

Here, we consider that the number of macule, papule, and vesicle lesions (i.e., L₁) increases at per lesion rate ge^−ct as the infection progresses, where g corresponds to the maximum growth rate of the number of lesions and c determines the rate at which this growth rate decreases over time due to viral clearance. It is assumed that the lesion counts in L₁ start to exponentially increase after t₀ days post first lesion onset and L₁ remains at 1 until that time (i.e., L₁(t) = 1 and L₂(t) = L₃(t) = L₄(t) = 0 for 0 < t ≤ t₀). The differential equations above therefore apply when the time since lesion onset t > t₀. In those equations, the number of scab lesions (i.e., L₄) decreases at per lesion rate d, and the transition rates between the stages are denoted by a_n(1 − e^−τ_nt) for n = 1,2,3 where a_n is the maximum transition rate and 1 − e^−τ_nt sets the time dependence of delays between the stage transitions. We estimated the above parameters, i.e., g, c, a_n, d, t₀, τ_n for n = 1,2,3, by model fitting (see below).

Modeling decay dynamics of viral load

To describe MPXV dynamics in peripheral blood, we used a simple decay model, which was previously used in a study of mpox [39]: where the variable V(t) is the viral RNA load (copies/mL) at time t and parameter δ represents the viral clearance rate. Note that the variable t again represents time after lesion onset; t = 0 is thus the date on which the first lesion was identified.

Parameter estimation

We estimated the parameters of both the lesion transition and viral load models using a nonlinear mixed-effect modeling approach, which incorporates fixed effects as well as random effects describing the inter-patient variability in parameter values. Each parameter of patient k, θ_k(= θ × e^πk) is decomposed as product of θ (the fixed effect) and e^πk (the random effect), where π_k is assumed to be drawn from a normal distribution: N(0, Ω). The fixed effect parameters and the standard deviations of random effects, Ω, were estimated using the stochastic approximation expectation/maximization (SAEM) algorithm, and individual parameters for each patient were then estimated using the empirical Bayes method. A right-truncated normal distribution was used in the likelihood function to account for the left censoring of the lesion count data or viral load data (i.e., when the lesion count or viral load is not detectable) [59]. MONOLIX 2021R2 (www.lixoft.com), a program for maximum likelihood estimation for nonlinear mixed-effect models, was used to fit the model to these data. We changed the initial values multiple times to avoid obtaining a local minimum of the negative log-likelihood function and to confirm the robustness of our parameter estimates.

Unsupervised clustering and stratification of lesion transition dynamics

We extracted three key “features” of lesion transition dynamics for each individual: i) the average macule, papule, and vesicle lesion count increase rate over 40 days (i.e., T = 40); ii) the average stage transition rates for n = 1,2,3), again with T = 40; and iii) the scab lesion decrease rate, d. For the purpose of stratification of the lesion transition dynamics, unsupervised random forest clustering based on each of five features was used (rfUtilities in R). The use of random forest allowed us to avoid overfitting by bootstrap aggregating (bagging) and to achieve better generalization performance [60]. Specifically, after a random forest dissimilarity (i.e., the distance matrix between all pairs of samples) was obtained, it was visualized by using UMAP in a two-dimensional plane and was stratified with spectral clustering (Python scikit-learn). The optimal number of clusters was determined using the eigengap heuristic method [61].

Classifying the stratified groups with a logistic regression

We built a predictive model for the stratified groups using a logistic regression. In this model, we used the viral load at lesion onset as the predictor variable, with the classification of the stratified group as the outcome variable. To ensure that the model was both robust and generalizable, we subjected it to rigorous validation through a 6-fold cross-validation process. We assessed the model performance using AUC scores and reported the mean and standard deviation of these scores obtained from 50 replications [46].

Classifying the stratified groups with a LightGBM model

Since some features showed high correlation or had missing values, we used the gradient boosting machine learning algorithm, specifically LightGBM, as our primary method of analysis [62]. LightGBM was chosen for its decision tree-based algorithm, which inherently handles the multicollinearity of predictors effectively. Given the limited number of data samples, we did not fine-tune the hyperparameters of LightGBM and the model was constructed using its default settings.

Statistical analysis

When necessary, variables were compared among different groups using analysis of variance (the Kruskal-Wallis test and Fisher’s exact test, for more than two numerical or categorial variables, respectively) with false discovery rate, or the Mann-Whitney U test (for two numerical variables) with Bonferroni correction. All statistical analyses were performed using R (version 4.2.0).