This work was conducted as part of the West African Centre for Cell Biology of Infectious Pathogens (WACCBIP) Long-term HIV Infection Cohort (WHICH Study) (9). A longitudinal design was employed for ART-naïve (M0) participants—PLWH who were yet to start ART and then followed up at six- and twelve-months post ART. In parallel, a cross-sectional design was used for ART-experienced participants (T_E0) who had received ART for at least six months at enrolment. Recruitment was conducted between July 2022 and September 2024 at Greater Accra Regional Hospital, University of Ghana Hospital–Legon, Tema General Hospital, Ho Municipal Hospital, Upper East Regional Hospital (Bolgatanga), and War Memorial Hospital (Navrongo). Healthy controls (CON) were recruited from the International Maritime Hospital and the West African Centre for Cell Biology of Infectious Pathogens. In total, the study enrolled 32 healthy controls, 141 ART-naïve, 52 at six months, 23 at twelve months, and 74 ART-experienced.

Venous blood (10 mL) was collected into BD Vacutainer^® K2EDTA and SST™ tubes (BD Biosciences, UK). Plasma and serum were separated by centrifugation (2500 rpm for 10 minutes) and stored at –80°C. Peripheral blood mononuclear cells (PBMCs) and red blood cells (RBCs) were also isolated and cryopreserved.

Plasma viral RNA was extracted using the Quick-RNA Viral Kit (Zymo Research, Cat. No. R1035) following the manufacturer’s protocol. HIV-1 viral load was quantified using the Bosphore^® HIV-1 Quantification Kit (Anatolia Geneworks, Cat. No. ABHIQ3) on the QuantStudio™ 5 Real-Time PCR System (Applied Biosystems). Viral load values, initially obtained in International Units/mL (IU/mL), were converted to copies/mL using a conversion factor of 1 IU = 0.7 copies/ml, as specified by the manufacturer, this was to ensure comparability with the WHO International Standard for HIV RNA NAT assays (NBSIC code 97/650).

Plasma concentration of various cytokines and chemokines were evaluated using the Human Cytokine Magnetic 25-Plex Panel (Invitrogen, Thermo Fisher Scientific, USA). The cytokines assessed were granulocyte-macrophage colony-stimulating factor (GM-CSF), interferon alpha (IFN-α), interferon beta (IFN-β), interleukin-1 receptor antagonist (IL-1Ra), IL-1 beta (IL-1β), IL-2, IL-2 receptor (IL-2R), IL-4, IL-5, IL-6, IL-8 (CXCL8), IL-10, IL-12/IL-23p40, IL-13, IL-15, IL-17A, and tumor necrosis factor alpha (TNF-α). The chemokines were regulated on activation, normal T cell expressed and secreted (RANTES/CCL5), macrophage inflammatory protein-1 alpha (MIP-1α/CCL3), Eotaxin (CCL11), macrophage inflammatory protein-1 beta (MIP-1β/CCL4), monocyte chemoattractant protein-1 (MCP-1/CCL2), monokine induced by gamma interferon (MIG/CXCL9) and interferon gamma-induced protein 10 (IP-10/CXCL10).

The assay was conducted following the manufacturer’s instructions and as previously reported by Tapela et al. (14). Briefly, in a 96-well plate, 25 µL of antibody-coated beads were added and washed. Then, 100 µL of samples, standards, and blanks were added and incubated for 2 hours with shaking at 250 rpm on a MicroPlate Shaker (Thermo Scientific, Korea). Subsequently, 100 µL of biotinylated detector antibody was added and incubated for 1 hour. Thereafter 100 µL of streptavidin-RPE was added incubated for 30 minutes, wells were washed, and 150 µL of wash buffer was added. The assay was read using a Luminex MAGPIX system (Luminex Corporation, Austin, TX, USA) and data analyzed using xPONENT™ software (v4.3.229), according to the manufacturer’s protocol.

Cox proportional hazards regression was used to evaluate associations between cytokine levels and HIV progression at baseline, six months, and twelve months. Hazard ratios (HRs) and 95% confidence intervals (CIs) were estimated with the survival package in R and visualized with forest plots generated using the survminer package (15). Model significance was assessed with the log-rank test.

To identify cytokine predictors of viral load, we applied three machine learning modelling approaches: Least Absolute Shrinkage and Selection Operator (LASSO) regression to select cytokines with strong linear associations, Random Forest (RF) to capture non-linear interactions and estimate variable importance based on percentage increase in mean squared error (%IncMSE) and Gradient Boosting Machine (GBM) to assess relative influence across sequential decision trees (16, 17). Model performance was evaluated by comparing predicted versus observed log viral load on the test set. Variable importance plots identified top predictors. Partial dependence plots (PDPs) were used to visualize non-linear effects of the top 10 cytokines.

Cytokine–protein interaction networks were constructed using the STRING database (STRINGdb) (18); https://string-db.org/, focusing on cytokines differentially expressed between ART-naïve and ART-experienced groups. To ensure high-confidence interactions, only experimentally validated and high-confidence STRING interactions (confidence score > 0.7) were used.

Network structure was analyzed by computing key centrality measures, including degree centrality and betweenness centrality, to identify highly interconnected and functionally influential cytokines. The network was visualized using the igraph package in R (19). Functional enrichment was performed with clusterProfiler (20), using Gene Ontology (GO) (biological process, molecular function, cellular component) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (21, 22). Gene names were retrieved from UniProt (23), and terms with false discovery rate (FDR) < 0.05 were considered significant.

Cytokines and chemokines were categorized into three functional groups based on their established roles in HIV pathogenesis. Pro-inflammatory cytokines were IL-1β, IL-5, IL-6, TNF-α, IL-12/IL-23p40, GM-CSF, IFN-γ, IFN-α, IL-2, IL-7, IL-15, IL-2R, and IL-17. Anti-inflammatory; IL-10, IL-1Ra, IL-4 and IL-13. Chemokines were MIP-1α (CCL3), MIP-1β (CCL4), RANTES, Eotaxin, MCP-1, IL-8, MIG, and IP-10 (24–26). Raw data was processed in Microsoft Excel and analyzed using GraphPad prism software Inc version 8 (GraphPad Software, San Diego, CA, USA) and open resource packages anchored in R software version 4.1.0 (R Development Core Team, Vienna, Austria, and R studio Version 2024.12.0.467). Cytokine and chemokine concentrations were expressed as Net Median Fluorescence Intensity (net MFI). Viral load, cytokine, and chemokine data were log₁₀-transformed. Group comparisons were made using the Kruskal–Wallis test with Dunn’s post hoc test. Spearman’s correlations assessed associations between cytokines and viral load. Significance was set at p < 0.05.

HIV-1 viral load remained high and unsuppressed among ART-naive their longitudinal follow-up pairs (Table 1). In contrast, viral suppression was observed in ART-experienced participants at the time of recruitment, the majority of whom had been on treatment for more than five years.

ART-naïve participants had significantly elevated levels of several pro-inflammatory cytokines, particularly IL-6 and IL-12/IL-23p40, whereas IL-1β was significantly reduced compared to the control group (Figure 1). Additionally, cytokines involved in T-cell homeostasis and activation, including IL-15, IL-2, IL-2R, and IL-7, were significantly higher in ART-naïve individuals than in uninfected controls (Figure 1). In contrast, levels of the anti-inflammatory cytokine IL-1Ra were significantly lower in ART-naïve individuals (Figure 2). Chemotactic cytokines such as MCP-1, IP-10, and MIG were also significantly elevated in ART-naïve individuals compared to controls (Figure 3).

Following ART initiation, cytokine and chemokine profiles shifted markedly. In ART-experienced individuals, there was a significant reduction in pro-inflammatory cytokines including GM-CSF, IL-12/IL-23p40, IL-6, IL-15, IL-17, and TNF-α compared to ART-naïve participants (Figure 1). Notably, anti-inflammatory cytokines IL-10 and IL-1Ra were also significantly reduced in ART-experienced individuals compare to ART-naïve (Figure 2). Among chemokines, Eotaxin and RANTES were significantly increased in ART-experienced whereas MCP-1, IL-8, MIG, and IP-10 were significantly decreased compared to ART-naïve (Figure 3).

A Cox proportional hazards regression was used to access the predictive capacity of cytokines andchemokines for HIV disease progression among ART-naive, six and twelve months follow up participants (Supplementary Figure 1). The model demonstrated strong predictive performance (AIC = 679.94, concordance index = 0.84) with a highly significant global log-rank p-value (p = 7.07 × 10^-5). Notably, lower levels of IL-1β (HR = 0.111, 95% CI: 0.0231–0.54, p = 0.006) and IL-10 (HR = 0.176, 95% CI: 0.0551–0.56, p = 0.003) were significantly associated with a higher likelihood of achieving virologic control. In contrast, elevated levels of GM-CSF (HR = 2.992, 95% CI: 1.118–8.01, p = 0.029) were associated with uncontrolled viraemia.

To evaluate cytokine predictors of HIV viral load, we applied LASSO regression, Random Forest (RF), and Gradient Boosting Machine (GBM) models. All three models demonstrated moderate predictive performance, with predicted versus observed log viral load showing good calibration (Supplementary Figure 2). RANTES and IP-10 were the strongest predictors, followed by IL-12/IL-23p40, Eotaxin, and IL-7 (Figure 4). Partial dependence analyses further highlighted non-linear cytokine–viral load relationships. RANTES exhibited an inverse association with viral load, while IP-10 and IL-12/IL-23p40 displayed positive effects, and Eotaxin and IL-7 showed threshold-dependent influences (Supplementary Figure 3).

The cytokine interaction network in HIV pathogenesis revealed a complex web of interactions between key pro-inflammatory and immunoregulatory cytokines (Figure 5). IL-1Ra emerged as a central regulatory node, exhibiting the highest betweenness centrality. Other highly connected cytokines included TNF-α, IL-6, IL-17, and IL-10. Chemokines such as RANTES, MIG, and IP-10 formed strongly interconnected nodes, with Eotaxin also showing notable centrality. IFN-γ, IL-2, IL-7, GM-CSF, IL-1β, and MIP-1α were also integrated into the network (Figure 5).

Spearman correlation analysis identified IP-10 as the strongest positive correlate of viral load, while RANTES showed a significant negative correlation (ρ ≈ –0.45) (Figure 6). Additional positive correlations with viral load were observed for IL-12/IL-23p40, MIG, MCP-1, IL-6, IL-2R, IL-2, and IFN-α.

Gene ontology (GO) enrichment analysis revealed distinct immune processes altered between ART-naïve and ART-experienced groups (Figures 7, 8; Supplementary Figure 4). Upregulated biological processes included chronic inflammatory response, eosinophil migration, chemokine-mediated signaling, granulocyte chemotaxis, and antimicrobial humoral responses. In contrast, pathways related to leukocyte activation, differentiation, and lymphocyte activation were downregulated.

At the molecular function level, chemokine activity, cytokine receptor binding, and G protein–coupled receptor binding were enriched, along with phospholipase activator activity. Downregulated functions included cytokine activity, growth factor receptor binding, and CXCR chemokine receptor interactions. In the cellular component category, downregulated pathways were mainly associated with the external plasma membrane and receptor complexes.

KEGG pathway enrichment analysis of cytokines and chemokines that differed significantly between ART-naïve and ART-experienced groups revealed enrichment of immune pathways, including viral protein–cytokine interactions, chemokine signaling, and cytokine–cytokine receptor interactions (Figure 9). Downregulated pathways included IL-17 and JAK-STAT signaling, along with those related to inflammatory bowel disease, rheumatoid arthritis, malaria, Chagas disease, hematopoietic lineage differentiation, and allograft rejection.