Machine learning aided multiscale modelling of the HIV-1 infection in the presence of NRTI therapy

An artificial neural network model for isolate-fold change relation

There exists various genotype-phenotype experiment data, including the fold change values of IC₅₀ (the required drug concentration to inhibit 50% of virions) for various reverse transcriptase inhibitors in the presence of susceptible and resistant isolates (Rhee et al., 2005). The most used genotype-phenotype data is the Stanford HIV drug resistance database (https://hivdb.stanford.edu/). We use filtered genotype-phenotype data of reverse transcriptase inhibitors available in this database and are widely used for various machine learning algorithms (Amamuddy, Bishop & Bishop, 2017; Masso & Vaisman, 2013). By regulating the data for each NRTI, 1,224 unique mutations were observed for the reverse transcriptase enzyme. In this filtered dataset, 1,662 isolates for epivir (3TC), 1,597 isolates for abacavir (ABC), 1,683 isolates for zidovudine (AZT), 1,693 isolates for stavudin (D4T), 1,693 isolates for didanosine (DDI) and 1,354 isolates for tenofovir (TDF) have been analyzed for NRTI susceptibility. The dataset includes 1,206, 1,136, 1,220, 1,223, 1,223, and 1,119 unique mutations for 3TC, ABC, AZT, D4T, DDI, and TDF, respectively.

Here, we apply the binary barcoding technique (Rhee, Taylor & Fessel, 2010) to represent the isolates occurring in the dataset. Hence, 1,224-dimensional input vectors of 0s and 1s are created by considering the existence of unique mutations in the isolates. Let us denote our complete mutation set as $M=\left\{m_1,m_2,\dots ,m_{1224}\right\}$ where m_i is an NRTI specified mutation pattern. We define the binary representation of isolate j as $I_j=\left[a_1,a_2,\dots ,a_{1224}\right]$ with ${\displaystyle a_k=\left\{\begin{array}{c}\hfill 1,if{m}_k\in I_j\hfill \\ \hfill 0,otherwise.\hfill \end{array}\right.}$ We construct six artificial neural networks (ANN) models to predict logarithmic fold-change values in the presence of any isolates related to each NRTI therapy by using the Machine Learning and Deep Learning toolbox of the MATLAB 2022a program (https://www.mathworks.com/). The ANN architectures include 1,224-dimensional input, five hidden layer neurons, and one output neuron with hyperbolic tangent-sigmoid and linear activation functions. The model selection process is explained with detailed quantitative observations in Table S1. The scaled conjugate gradient algorithm with MATLAB built-in function “trainscg” has been used in the training process over GPU. Let us denote our model as a function that maps isolate vectors to the fold changes as ${\displaystyle FoldChange={ANN}_X\left(I\right)}$ where I ∈ {0, 1}^1×1224 and X is a specified inhibitor (X ∈ {3TC, ABC, AZT, D4T, DDI, TDF}). To overcome possible overfitting, we have implemented an ensemble learning process. For each inhibitor, the 50 ×100 model has been trained with random training, validation, and test set (80%, 10% and 10%). A model is chosen from every 100 models that yield the minimum mean square error for the test set of the corresponding inhibitor data. Hence, 50 optimal models are selected out of 5,000 models for each NRTI inhibitor, and the final model is calculated as the average of these models.

The prediction performance of six ANN_X(I) models with linear correlation coefficient (R) and mean square error (MSE) values are presented in Fig. 1. According to the figure, ANN_X(I) models yield accurate predictions with high R and low MSE scores. Mean MSE value of ANN_X(I) models have been obtained as 0.0453 with 95% CI [0.0005–0.0901]. Similarly, the mean R value of the models has been calculated as 0.9093 with 95% CI [0.8677–0.9509]. To observe how six ANN_X(I) models classify resistant and susceptible strains, we convert our regression models into classification models by labeling the data as resistant (Fold Change ≥ 3) and susceptible (Fold Change < 3). The receiving operating curves (ROC) corresponding to the six ANN models and the area under the curve (AUC) values are presented in Fig. S1. According to the classification results, we get the mean AUC score as 0.9649 with 95% CI [0.9423–0.9875]. Additionally, to see why such a nonlinear model is needed to map the genotype data into the phenotype output, we also perform multiple linear regression (MLR) analysis (with 20% holdout data) for data of six NRTIs. The regression and classification performance of the MLR models are shown in Figs. S2–S3. A fair comparison between the ANN and MLR models in terms of the MSE, R, and AUC values is given in Table S2. According to the table, even classification performance of the models is almost the same, the ANN models give much more accurate estimations in regression. Since better regression performance is more desirable for our further modelling framework, the ANN models are assumed to be our baseline models for predicting the resistance profiles of given viral strains.

Modelling the time-dependent drug efficacy

Modelling the efficacy of antiretrovirals using the plasma drug concentrations can be seen in various studies in the literature (Dixit & Perelson, 2004; Rong, Feng & Perelson, 2007; Rosenbloom et al., 2012). Rosenbloom et al. (2012) modeled the time-dependent drug efficiency in plasma by considering the exponential decay of plasma drug concentration after the instantaneous peak. Here we use the time-dependent drug efficacy model described by Dixit & Perelson (2004) and Rong, Feng & Perelson (2007), considering the pharmacokinetic parameters of drugs in the blood. Dixit & Perelson (2004) considered the phosphorylated concentration of the tenofovir (TDF) in the cells. Since the time-drug efficiency functions obtained by taking into account blood concentration and phosphorylated within cell concentration of drugs follow a very similar trend, here we assume the blood concentration of the drugs (see Fig. 1 of Dixit & Perelson (2004)). Additionally, the non-availability of phosphorylation reaction parameters for the remaining five inhibitors 3TC, ABC, AZT, D4T, and DDI have encouraged us to consider the blood concentration of the drugs only.

Let ${\varepsilon }_X^Y\left(t\right)$ denotes the time-dependent efficacy of drug X in the presence of strain (isolate) Y. The instantaneous efficacy can be approximated as Dixit & Perelson (2004) (2)${\displaystyle {\varepsilon }_X^Y\left(t\right)=\frac{C_b^X\left(t\right)}{{\left({IC}_{50}\right)}_X^Y+C_b^X\left(t\right)}}$ where $C_b^X\left(t\right)$ denotes the within blood concentration of drug X and ${\left({IC}_{50}\right)}_X^Y$ denotes the required concentration of drug X to inhibit the 50% of strain Y. According to our isolate-fold change ANN model, Eq. (2) can be rewritten as (3)${\displaystyle {\varepsilon }_X^Y\left(t\right)=\frac{C_b^X\left(t\right)}{{ANN}_X\left(Y\right){\left({IC}_{50}\right)}_X^{WT}+C_b^X\left(t\right)}}$ where ${\left({IC}_{50}\right)}_X^{WT}$ denotes the required concentration of drug X to inhibit the 50% wild type virus. Thus, to completely describe ${\varepsilon }_X^Y\left(t\right),$ we should model $C_b^X\left(t\right).$ According to Dixit & Perelson (2004), the concentration of a drug in the blood can be expressed as (4)${\displaystyle C_b\left(t\right)=\frac{FD{k}_a{e}^{-k_{e}t}}{V_d\left(k_e-k_a\right)\left(e^{k_a{I}_d}-1\right)}\left[\right.1-e^{\left(k_e-k_a\right)t}\left(1{-e}^{N_d{k}_a{I}_d}\right)+\frac{\left(e^{k_e{I}_d}-e^{k_a{I}_d}\right)\left(e^{{\left(N_d-1\right)k}_e{I}_d}-1\right)}{e^{k_e{I}_d}-1}-e^{\left({\left(N_d-1\right)k}_e+k_a\right)I_d}\left]\right.}$

where F is the bioavailability of the drug, D is the mass of the drug administered in one dose, I_d is the dosing interval, N_d is the number of doses up to time t, V_d is the volume of distribution, k_a and k_e are pharmacokinetic parameters. The drug-specific parameters k_a,  k_e,  D,  I_d and F occurred in Eq. (4) and IC₅₀ values for 3TC, ABC, AZT, D4T, DDI, and TDF according to the equations given by Dixit & Perelson (2004) are evaluated and presented in Table 1. Detailed explanations of the derivation of these parameters are given in the Supplementary Information.

A multi-strain within-host model

This part of the study combines our investigations into a unique multi-strain within-host model. To reduce the cost of the simulations, we assume the main NRTI-related mutations 115F, 151M, 184I, 184V, 210W, 215F, 215Y, 41L, 65N, 65R, 67N, 69D, 70E, 70G, 70R, 74I and 74V according to the study of Rhee et al. (2005). These 17 mutations yield 131,071 unique strains having all possible mutations. Thus, by considering wild-type and mutant strains, we have total N = 131,072 strains. Our multi-strain within-host model with time-dependent NRTI therapy can be derived from one strain model (1) as follows ${\displaystyle \frac{dT}{dt}=s_T-k_{T}T\sum _{i=1}^N\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)V_i-{\delta }_{T}T+\frac{{\rho }_T\sum _{i=1}^N{V}_i}{c_T+\sum _{i=1}^N{V}_i}T}$ ${\displaystyle \frac{d{T}_i^{\ast }}{dt}=k_T\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)T{V}_i-{\delta }_{T^{\ast }}{T}_i^{\ast }}$ (5)${\displaystyle \frac{dM}{dt}=s_M-k_{MM}\sum _{i=1}^N\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)V_i-{\delta }_{M}M+\frac{{\rho }_M\sum _{i=1}^N{V}_i}{c_M+\sum _{i=1}^N{V}_i}M}$ ${\displaystyle \frac{d{M}_i^{\ast }}{dt}=k_M\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)M{V}_i-{\delta }_{M^{\ast }}{M}^{\ast }}$ ${\displaystyle \frac{d{V}_i}{dt}=p_T{T}_i^{\ast }+p_M{M}_i^{\ast }-{\delta }_V{V}_i}$

where i = 1, 2, …, N =131,072, $T_i^{\ast }\left(t\right)$ and $M_i^{\ast }\left(t\right)$ denote the number of CD4 + T cells and macrophage cells infected by strain i and V_i(t) represents the number of virions having i th genotype. In the multi-strain within-host model (5), ${\varepsilon }_X^i\left(t\right)$ denotes the time-dependent efficacy of the inhibitor X on the strain i and 0 ≤ c_i ≤ 1 represents the fitness costs of mutant strains with c₁ = 1 for the wild type of strain. The lack of enough experimental results on these fitness values compelled us to use the mean fitness cost values of mutations 41L, 67N, 70R, 184V, 210W, 215D, 215S, and 219Q estimated by Kühnert et al. (2018) as 0.2232, 0.3181, 0.3863, 0.5899, 0.3091, 0.0981, 0.1664 and 0.3207, respectively. According to these data, we assume that c_i = 0.3015 for mutant strains i ≥ 2. A schematic illustration of the multi-strain within-host model (5) is given in Fig. 2. Parameter values of multi-strain within-host model (5) with corresponding references can be seen in Table 2.

The within-host model (5) ignores the role of latently infected CD4+ T cells. The main role of latently infected CD4+ T cells is the viral rebound after poor adherence to the given therapy (Chun et al., 2000), and these cells are almost three percent of all CD4+ T cells (Hadjiandreou, Conejeros & Wilson, 2009). Since model (5) is continuous over time and hence the emerged viral strains are not completely eradicated in the viral suppression phase, the persistence of HIV-1 virions is automatically ensured, and poor adherence in model (5) provides viral rebound. Thus, ignoring the latently infected CD4+ T cells in model (5) does not considerably affect our modelling framework. As indicated in the study of Chun et al. (2000), latently infected CD4+ T cells are not the only reason for the rebound of plasma viremia after discontinuation of the ART. The literature (Alexaki, Liu & Wigdahl, 2008; Hendricks et al., 2021; Kruize & Kootstra, 2019) show that the macrophage cells are of particular importance in HIV-1 persistence, and this is why model (5) considers this observation like some existing studies (Hadjiandreou, Conejeros & Vassiliadis, 2007; Hadjiandreou, Conejeros & Wilson, 2009; Hernandez-Vargas, 2019; Hernandez-Vargas & Middleton, 2013). Consideration of the role of the macrophage cells yields slow progression to the AIDS phase for untreated patients and improves the reliability of model outcomes (Hernandez-Vargas, 2019).

To model the effect of mutations, we do not explicitly include the mutation matrix in the ODE system (5); instead, we address the transition between mutations and strains at the end of each time step by generating Poisson random numbers (Rosenbloom et al., 2012). Let us assume time step n (t = n day), ${\left(T_i^{\ast }\right)}_n=T_i^{\ast }\left(n\right)$ and ${\left(M_i^{\ast }\right)}_n=M_i^{\ast }\left(n\right)$. The mutation matrix of our system is denoted by MT and defined as (6)${\displaystyle {MT}_{ij}=\left\{\begin{array}{c}\hfill 1,ifstrainicantakeamutationtobecomestrainj\hfill \\ \hfill 0,otherwise\hfill \end{array}\right.}$

For the infected CD4 T cells ${\left(T_i^{\ast }\right)}_n$ and infected macrophage cells ${\left(M_i^{\ast }\right)}_n$, we calculate the number of new infected ones in one day period as $\Delta {\left(T_i^{\ast }\right)}_n$ and $\Delta {\left(M_i^{\ast }\right)}_n$ without taking into account the death of these newly infected cells. For each i = 1, 2, …, N, $poissrnd\left(\mu \Delta {\left(T_i^{\ast }\right)}_n\right)$ and $poissrnd\left(\mu \Delta {\left(T_i^{\ast }\right)}_n\right)$ number of infected cells are randomly transmitted from strain i to strain j according to the mutation matrix MT_ij where function poissrnd(x) generates Poisson random number with mean x and μ = 3 × 10⁻⁵ denoting the mutation rate (Rosenbloom et al., 2012). Note that the mutation rate for each point mutation is unique for the corresponding amino acid change, but we assume a fixed average mutation rate μ = 3 × 10⁻⁵ as stated by Rosenbloom et al. (2012). Since NRTI-related mutation rates have low variance value (Rosenbloom et al., 2012) and we have so many viral strains to track, we use overall mutation rate μ = 3 × 10⁻⁵. Parameter values of models (1) and (5) are presented with their references in Table 2.

Model (5) can also include dual therapy of NRTIs X and Y by modifying the therapy-related time-dependent infection coefficients for CD4 + T cells and macrophage cells ${\beta }_i^{T/M}\left(t\right)=k_{T/M}\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)$ with the use of Bliss independence of drug actions as Jilek et al. (2012) (7)${\displaystyle {\beta }_i^{T/M}\left({\varepsilon }_X^i\left(t\right),{\varepsilon }_Y^i\left(t\right)\right)=k_{T/M}\left(1-c_i\right)\left(1-{\varepsilon }_X^i\left(t\right)\right)\left(1-{\varepsilon }_Y^i\left(t\right)\right)}$ or Loewe additivity of drug actions (Jilek et al., 2012) (8)${\displaystyle {\beta }_i^{T/M}\left({\varepsilon }_X^i\left(t\right),{\varepsilon }_Y^i\left(t\right)\right)=k_{T/M}\left(1-c_i\right)\frac{1}{\frac{{\varepsilon }_X^i\left(t\right)}{1-{\varepsilon }_X^i\left(t\right)}+\frac{{\varepsilon }_Y^i\left(t\right)}{1-{\varepsilon }_Y^i\left(t\right)}+1}.}$

Bliss independence assumes independent actions of combined drugs, and Loewe additivity assumes the competition for the same binding site. According to Jilek et al. (2012), all combinations except AZT-D4T and DDI-TDF obey the Bliss independence rule, and these two combinations obey the Loewe additivity rule. Note that, since we assume k_M ≈ k_T/1000 and ${\beta }_i^T\left(t\right)\approx {\beta }_i^M\left(t\right)/1000$ according to the Hernandez-Vargas (2019) and Hernandez-Vargas & Middleton (2013) (see Table 2), we prefer to use the notation β_i for ${\beta }_i^T$ throughout the following parts. Whenever β_i values are quantitatively mentioned in the results section, these values correspond to the ${\beta }_i^T.$

Note that even though we describe our model parameters for 1 ml of blood in Table 2 as widely assumed in the literature (Hadjiandreou, Conejeros & Vassiliadis, 2007; Hernandez-Vargas, 2019; Hernandez-Vargas & Middleton, 2013), we simulate the viral dynamics in the host plasma (3,000 ml; Rosenbloom et al., 2012) to catch more viral diversity. We assume that the only reservoir of HIV virions is the plasma, which is the major one (Valcour et al., 2012), even if there exist other reservoirs like lymph nodes or cerebrospinal fluid (CSF) (Valcour et al., 2012; Haase, 1999). Since the instantaneous drug efficiency rates are (${\varepsilon }_X^Y\left(t\right)$) in non-dimensionless form, we can easily simulate the dynamics in the host plasma by converting the volume-dependent model parameters given in Table 2. For example, by considering 3L host plasma (Rosenbloom et al., 2012), the infectivity parameter k_T = 4.5714 × 10⁻⁸ ml/day equivalently becomes $k_T=\frac{4.5714\times {10}^{-8}}{3000}$ plasma/day = 1.5238 × 10⁻¹¹ plasma/day.