Machine learning aided multiscale modelling of the HIV1 infection in the presence of NRTI therapy
 Published
 Accepted
 Received
 Academic Editor
 Yuriy Orlov
 Subject Areas
 Computational Biology, Mathematical Biology, HIV, Statistics
 Keywords
 AIDS, HIV infection, Machine learning, NRTI therapy, Mathematical models
 Copyright
 © 2023 Tunc et al.
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
 Cite this article
 2023. Machine learning aided multiscale modelling of the HIV1 infection in the presence of NRTI therapy. PeerJ 11:e15033 https://doi.org/10.7717/peerj.15033
Abstract
Human Immunodeficiency Virus (HIV) is one of the most common chronic infectious diseases in humans. Extending the expected lifetime of patients depends on the use of optimal antiretroviral therapies. Emergence of the drugresistant strains can reduce the effectiveness of treatments and lead to Acquired Immunodeficiency Syndrome (AIDS), even with antiretroviral therapy. Investigating the genotypephenotype relationship is a crucial process for optimizing the therapy protocols of the patients. Here, a mathematical modelling framework is proposed to address the impact of existing mutations, timing of initiation, and adherence levels of nucleotide reverse transcriptase inhibitors (NRTIs) on the evolutionary dynamics of the virus strains. For the first time, the existing Stanford HIV drug resistance data have been combined with a multistrain withinhost ordinary differential equation (ODE) model to track the dynamics of the most common NRTIresistant strains. Overall, the D4T3TC, D4TAZT and TDFD4T drug combinations have been shown to provide higher success rates in preventing treatment failure and further drug resistance. The results are in line with the genotypephenotype data and pharmacokinetic parameters of the NRTI inhibitors. Moreover, we show that the undetectable mutant strains at the diagnosis have a significant effect on the success/failure rates of the NRTI treatments. Predictions on undetectable strains through our multistrain withinhost model yielded the possible role of viral evolution on the treatment outcomes. It has been recognized that the improvement of multiscale models can contribute to the understanding of the evolutionary dynamics, and treatment options, and potentially increase the reliability of genotypephenotype models.
Introduction
Antiretroviral drug resistance is one of the main barriers to therapy success for HIVpositive patients. According to the WHO, the HIV drug resistance report 2021, 10% and 40% of adults are affected by drugresistant strains (DRS) for naive and treated patients, respectively. In addition, 50% of newly diagnosed infants were exposed to the DRS. The DRS can be acquired with nonadherence to the therapy protocols, or patients can directly be infected with DRSs (Blower et al., 2001). Both scenarios yield lifelong persistence of the DRS and need to be carefully tracked.
Quantitative evaluation of HIV drug resistance has been carried out with the use of phenosense assays by finding the foldchange of IC_{50} values (the amount of concentration to inhibit 50% of virion) between drugresistant and wildtype strains (Zhang et al., 2005; Pham et al., 2018; Feng et al., 2016). Data modelling frameworks have been used to construct general mathematical relations between genotype and phenotype information (Tarasova et al., 2018; Steiner, Gibson & Crandall, 2020; Shah et al., 2020; Tarasova et al., 2023; Lagunin et al., 2023). These mathematical models aim to generalize the given data by means of encoding the amino acid sequence of target enzymes (Rhee, Taylor & Fessel, 2010). One of the main contributions of the current study is to explore how these models can be embedded into a withinhost model to simulate the evolutionary dynamics of HIV strains. In particular, we simulate thousands of clinically relevant HIV mutant strains provided by Rhee, Taylor & Fessel (2010).
For forecasting the viral dynamics of HIV, various withinhost models have been presented in ordinary differential equation (ODE) forms in the presence/absence of resistant strains and antiretroviral therapy (Perelson & Nelson, 1999; Dixit & Perelson, 2004; Rong, Feng & Perelson, 2007; Hadjiandreou, Conejeros & Vassiliadis, 2007; Sutimin et al., 2017; Wu & Zhao, 2020; Chen, Teng & Zhang, 2021). Perelson & Nelson (1999) proposed HIV withinhost models consisting of CD4+ T cells (T), infected CD4+ T cells (T^{∗}), macrophage cells (M), infected macrophage cells (M^{∗}) and virions (V) in the presence and absence of antiretroviral therapy. Dixit & Perelson (2004) proposed T—T^{∗}—V model by considering timedependent intracellular efficiency of reverse transcriptase and protease inhibitors (RTIs and PIs). Rong, Feng & Perelson (2007) derived twostrain extension of the withinhost model given in the literature (Perelson & Nelson, 1999; Dixit & Perelson, 2004) with antiretroviral therapy. Hadjiandreou, Conejeros & Vassiliadis (2007) revised the T—T^{∗}—M—M^{∗}—V model proposed by Perelson & Nelson (1999) by adding homeostatic cell proliferation terms to capture long time behaviour of the HIV dynamics. Sutimin et al. (2017) modeled the withinhost HIV dynamics with target Langerhans and CD4+ T cells and investigated the timedependent efficiency of RTIs and PIs with various scenarios. Wu & Zhao (2020) derived twostrain withinhost model including the age of infection detail represented by the system of integrodifferential equations. They mathematically formulated the competition between the drugsensitive and drugresistance strains with respect to model parameters. Chen, Teng & Zhang (2021) included saturated incidence and distributed infection delays into the standard twostrain T—T^{∗}—V model and investigated the effects of those novel incidence terms on the longtime behaviour of the dynamics. Additionally, the effect of drug adherence on the virological failure of ARTs (Rosenbloom et al., 2012), the effect of timedependent drug efficiencies on ART response (Rong, Feng & Perelson, 2007; Vaidya & Rong, 2017), competition between susceptible and resistant strains in the viral dynamics (Ball, Gilchrist & Coombs, 2007; Lythgoe, Pellis & Fraser, 2013), the role of latently infected CD4+ T cell reservoirs on the evolution of strains (Doekes, Fraser & Lythgoe, 2017), comparison of entry inhibitors with the RTIs and PIs according to viral resistance (Alshorman, Alhosainat & Jackson, 2022), investigation of optimal timing for ART Rouzine (2022) have been proposed through withinhost models. The proposed mathematical models assume the coexistence of susceptible and resistant strains and generally investigate the response to antiretroviral therapy (ART). The current study addresses similar questions with a novel multiscale model based on Stanford HIV Drug Resistance data and machine learning models.
For the first time, we combined the experimental drug resistance data of nucleotidereverse transcriptase inhibitors (NRTI) available in the Stanford HIV drug resistance database (Rhee et al., 2003) with a withinhost model of HIV infection to observe the dynamics of the viral strains under different scenarios. Our multiscale model brings together three pieces of information: IC_{50} values for each mutant with machine learning models, within blood dynamics for NRTIs, and CD4+ T cells and macrophage cells for primary targets of virions. For different mutant compositions, we aim to investigate the emergence of treatment failure for different initiation timing (up to one year) and adherence level of NRTI therapies (21 different combinations). Here we rank the inhibitory capabilities of the NRTI combinations in the presence of various viral strains and ongoing viral evolution. Our results add to the predictions of the Stanford HIV drug resistance database, which identifies the best drug by selecting the one that has the lowest IC_{50} for a given mutant. But that model is a static model that cannot incorporate the effects of new mutants that can be generated through time which is accounted for in our model.
Materials and Methods
Withinhost model with wildtype virus
In this part, we have inspired from the earlier studies on the withinhost HIV infection model (Hadjiandreou, Conejeros & Vassiliadis, 2007; HernandezVargas, 2019; HernandezVargas & Middleton, 2013). We assume that the primary reservoirs for HIV infection are: CD4+T cells and macrophages denoted by T(t) and $M\left(t\right)$ (HernandezVargas, 2019; HernandezVargas & Middleton, 2013). The longliving macrophage cells cause the persistence of virions over the years (Orenstein, 2001; Herbein & Varin, 2010). Macrophage cells contribute to the depletion of healthy CD4 + T cells in advanced HIV infection (Crowe, 1995). Withinhost modelling of HIV infection without considering the macrophage reservoirs yielded less reliable dynamics, such as the models that never result in the AIDS phase (Rong, Feng & Perelson, 2007). We denote the HIV infected CD4+ T cells and macrophages by T^{∗}(t) and ${M}^{\ast}\left(t\right).$ Lastly, the number of free wildtype virions in the host is denoted by the function $V\left(t\right).$ By considering model assumptions like homeostatic cell proliferation terms (s_{T}, s_{M}), bilinear incidence terms (k_{T}TV, k_{M}TM), natural deaths of cells and virions (δ_{T}T, δ_{M}M, δ_{T∗}T^{∗}, δ_{M∗}M^{∗}, δ_{V}V), viral replication terms (p_{T}T^{∗}, p_{M}M^{∗}) and the Michaelis–Menten type proliferation terms $\left(\frac{{\rho}_{T}V}{{c}_{T}+V}T,\frac{{\rho}_{M}V}{{c}_{M}+V}M\phantom{\rule{0.25em}{0ex}}\right)$, we express the one strain withinhost model with the following system of ordinary differential equations (HernandezVargas, 2019; HernandezVargas & Middleton, 2013) $\frac{dT}{dt}={s}_{T}{k}_{T}TV{\delta}_{T}T+\frac{{\rho}_{T}V}{{c}_{T}+V}T$ $\frac{d{T}^{\ast}}{dt}={k}_{T}TV{\delta}_{{T}^{\ast}}{T}^{\ast}$ (1)$\frac{dM}{dt}={s}_{M}{k}_{M}MV{\delta}_{M}M+\frac{{\rho}_{M}V}{{c}_{M}+V}M$ $\frac{d{M}^{\ast}}{dt}={k}_{M}MV{\delta}_{{M}^{\ast}}{M}^{\ast}$ $\frac{dV}{dt}={p}_{T}{T}^{\ast}+{p}_{M}{M}^{\ast}{\delta}_{V}V$ where initial conditions are considered as $T\left(0\right)={T}_{0},\phantom{\rule{0.25em}{0ex}}{T}^{\ast}\left(0\right)={T}_{0}^{\ast},\phantom{\rule{0.25em}{0ex}}M\left(0\right)={M}_{0},\phantom{\rule{0.25em}{0ex}}{M}^{\ast}\left(0\right)={M}_{0}^{\ast}$ and $V\left(0\right)={V}_{0}$. Further details of the model (1) can be seen in the study of HernandezVargas & Middleton (2013). In the following section, we expand the model Eq. (1) to include both susceptible and resistant multiple strains as well as NRTI therapy.
Multi strain withinhost model with NRTI therapy
The ARTs include at least one of the NRTIs that aim to block the activation of the reverse transcriptase enzyme. Effective treatment of HIVpositive patients with NRTIs saves millions of lives worldwide (Tressler & Godfrey, 2012). However, the errorprone structure of the HIV replication yields resistant strains over the years, and these strains are known to be a primary barrier to preventing AIDS (Kuritzkes, 2011). Our multiscale withinhost model includes three main steps: constructing machine learning models to generalize isolatefold change data for NRTIs, a model for dealing with NRTI action in blood, and finally, a withinhost model with multistrains and NRTI therapy.
An artificial neural network model for isolatefold change relation
There exists various genotypephenotype experiment data, including the fold change values of IC_{50} (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 genotypephenotype data is the Stanford HIV drug resistance database (https://hivdb.stanford.edu/). We use filtered genotypephenotype 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,224dimensional 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},\phantom{\rule{0.25em}{0ex}}{m}_{2},\phantom{\rule{0.25em}{0ex}}\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},\phantom{\rule{0.25em}{0ex}}{a}_{2},\phantom{\rule{0.25em}{0ex}}\dots ,{a}_{1224}\right]$ with ${a}_{k}=\left\{\begin{array}{c}\hfill 1,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}if\phantom{\rule{0.25em}{0ex}}{m}_{k}\in {I}_{j}\hfill \\ \hfill 0,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}otherwise.\hfill \end{array}\right.$ We construct six artificial neural networks (ANN) models to predict logarithmic foldchange 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,224dimensional input, five hidden layer neurons, and one output neuron with hyperbolic tangentsigmoid and linear activation functions. The model selection process is explained with detailed quantitative observations in Table S1. The scaled conjugate gradient algorithm with MATLAB builtin 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 $Fold\phantom{\rule{0.25em}{0ex}}Change={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 timedependent 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 timedependent drug efficiency in plasma by considering the exponential decay of plasma drug concentration after the instantaneous peak. Here we use the timedependent 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 timedrug 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 nonavailability 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 timedependent efficacy of drug X in the presence of strain (isolate) Y. The instantaneous efficacy can be approximated as Dixit & Perelson (2004) (2)$\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 isolatefold change ANN model, Eq. (2) can be rewritten as (3)$\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)${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 drugspecific parameters k_{a}, k_{e}, D, I_{d} and F occurred in Eq. (4) and IC_{50} 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 multistrain withinhost model
This part of the study combines our investigations into a unique multistrain withinhost model. To reduce the cost of the simulations, we assume the main NRTIrelated 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 wildtype and mutant strains, we have total N = 131,072 strains. Our multistrain withinhost model with timedependent NRTI therapy can be derived from one strain model (1) as follows $\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$ $\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)$\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$ $\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}$ $\frac{d{V}_{i}}{dt}={p}_{T}{T}_{i}^{\ast}+{p}_{M}{M}_{i}^{\ast}{\delta}_{V}{V}_{i}$
Parameter/Drug  3TC  ABC  AZT  D4T  DDI  TDF 

IC_{50}(×10^{−5} mg/ml)  3.97  132.64  1.87  4.25  113.11  16.24 
D (mg)  300  300  300  40  400  300 
I_{d} (day)  1  0.5  0.5  0.5  1  1 
F  0.86  0.83  0.64  0.86  0.42  0.39 
k_{a}  27.98  51.07  37.42  54.29  32.34  8.36 
k_{e}  3.44  8.52  14.25  7.84  47.30  16.58 
V_{d} (ml)  91,000  60,200  112,000  46,000  54,000  87,500 
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 multistrain withinhost model (5), ${\varepsilon}_{X}^{i}\left(t\right)$ denotes the timedependent efficacy of the inhibitor X on the strain i and 0 ≤ c_{i} ≤ 1 represents the fitness costs of mutant strains with c_{1} = 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 multistrain withinhost model (5) is given in Fig. 2. Parameter values of multistrain withinhost model (5) with corresponding references can be seen in Table 2.
The withinhost 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 HIV1 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 HIV1 persistence, and this is why model (5) considers this observation like some existing studies (Hadjiandreou, Conejeros & Vassiliadis, 2007; Hadjiandreou, Conejeros & Wilson, 2009; HernandezVargas, 2019; HernandezVargas & 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 (HernandezVargas, 2019).
Parameter  Value  Unit  Parameter variation 

s_{T}  10^{4}  ml^{−1}d^{−1}  7 × 10^{3} − 2 × 10^{4} 
s_{M}  150  ml^{−1}d^{−1}  100 − 300 
k_{T}  4.5714 × 10^{−8}  mld^{−1}  3.2 × 10^{−8} − 10^{−7} 
k_{M}  4.3333 × 10^{−11}  mld^{−1}  1.73 × 10^{−11} − 1.3 × 10^{−9} 
p_{T}  38  d^{−1}  30.4 − 114 
p_{M}  35  d^{−1}  22 − 132 
δ_{T}  0.01  d^{−1}  0.001 − 0.017 
δ_{T∗}  0.4  d^{−1}  0.1 − 0.45 
δ_{M}  0.001  d^{−1}  10^{−4} − 1.4 × 10^{−3} 
δ_{M∗}  0.001  d^{−1}  10^{−4} − 1.2 × 10^{−3} 
δ_{V}  2.4  d^{−1}  0.96 − 2.64 
ρ_{T}  0.01  d^{−1}  – 
ρ_{M}  0.003  d^{−1}  – 
c_{T}  3 × 10^{5}  ml^{−1}  – 
c_{M}  2.2 × 10^{5}  ml^{−1}   
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)${MT}_{ij}=\left\{\begin{array}{c}\hfill 1,\phantom{\rule{0.25em}{0ex}}if\phantom{\rule{0.25em}{0ex}}strain\phantom{\rule{0.25em}{0ex}}i\phantom{\rule{0.25em}{0ex}}can\phantom{\rule{0.25em}{0ex}}take\phantom{\rule{0.25em}{0ex}}a\phantom{\rule{0.25em}{0ex}}mutation\phantom{\rule{0.25em}{0ex}}to\phantom{\rule{0.25em}{0ex}}become\phantom{\rule{0.25em}{0ex}}strain\phantom{\rule{0.25em}{0ex}}j\hfill \\ \hfill 0,\phantom{\rule{0.25em}{0ex}}otherwise\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\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^{−5} 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^{−5} as stated by Rosenbloom et al. (2012). Since NRTIrelated 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^{−5}. 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 therapyrelated timedependent 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)${\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)${\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 AZTD4T and DDITDF 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 HernandezVargas (2019) and HernandezVargas & 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; HernandezVargas, 2019; HernandezVargas & 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 nondimensionless form, we can easily simulate the dynamics in the host plasma by converting the volumedependent 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^{−8} ml/day equivalently becomes ${k}_{T}=\frac{4.5714\times {10}^{8}}{3000}$ plasma/day = 1.5238 × 10^{−11} plasma/day.
Results
This section provides the simulation results of the multistrain withinhost model (5), starting with various viral strains. The effects of adherence levels and initiation timing of NRTI therapies on the progression of viral dynamics are investigated. This section includes four subsections in which we propose the statistics of the infection rates, details of model simulations, the quantitative measure for the therapy success, and the simulation results for various cases.
Statistics of infection rates
Before running the simulations to observe the failure/success distribution of each NRTI combination, we may predict the best possible therapy protocol through our pretrained machine learning model and the pharmacokinetic properties of the inhibitors. Obviously, as we infer from our model (5) and drugspecific timedependent infection rate ${\beta}_{i}\left({\varepsilon}_{X}^{i}\left(t\right),{\varepsilon}_{Y}^{i}\left(t\right)\right)$ (7)–(8), each viral strain has its infection rate and aims to be dominant by infecting more healthy cells. Since evaluation of ${\beta}_{i}\left({\varepsilon}_{X}^{i}\left(t\right),{\varepsilon}_{Y}^{i}\left(t\right)\right)$ is straightforward through Eqs. (7)–(8) and (3), we may have some prior estimates for the selection of the best therapy protocol. Distribution of 131,071 ${\beta}_{i}\left({\varepsilon}_{X}^{i},{\varepsilon}_{Y}^{i}\right)={\int}_{0}^{1}{\beta}_{i}\left({\varepsilon}_{X}^{i}\left(t\right),{\varepsilon}_{Y}^{i}\left(t\right)\right)dt$ values in the presence of 21 different mono and dual NRTI therapies are illustrated in Fig. 3. Descriptive statistic values of ${\beta}_{i}\left({\varepsilon}_{X}^{i},{\varepsilon}_{Y}^{i}\right)$ values for all combinations are presented in Table 3.
Drugs  Mean  Min  Max  Std  Median  Mode  Q_{1}  Q_{3} 

D4T3TC  1.290  0.160  2.069  0.363  1.295  0.16  1.029  1.576 
D4TAZT  1.370  0.427  2.358  0.442  1.368  0.427  1.021  1.722 
TDFD4T  1.403  0.604  2.373  0.371  1.382  0.604  1.109  1.679 
D4T  1.442  0.697  2.319  0.351  1.426  0.697  1.157  1.72 
AZT3TC  1.473  0.154  2.212  0.462  1.531  0.154  1.109  1.877 
D4TABC  1.525  0.695  2.405  0.374  1.514  0.695  1.221  1.826 
DDID4T  1.592  0.765  2.466  0.382  1.581  0.765  1.279  1.903 
TDFAZT  1.627  0.474  2.523  0.506  1.683  0.474  1.226  2.056 
AZTABC  1.755  0.551  2.529  0.503  1.844  0.551  1.363  2.194 
AZT  1.775  0.554  2.564  0.513  1.858  0.554  1.373  2.223 
DDIAZT  1.834  0.562  2.602  0.533  1.933  0.562  1.412  2.307 
TDF3TC  1.884  0.3  2.265  0.307  1.952  0.3  1.835  2.065 
3TC  1.965  0.29  2.173  0.339  2.114  0.29  2.000  2.133 
ABC3TC  2.030  0.274  2.318  0.359  2.173  0.274  2.025  2.225 
DDI3TC  2.155  0.323  2.373  0.356  2.305  0.323  2.201  2.325 
TDFABC  2.172  1.508  2.588  0.181  2.176  1.508  2.038  2.314 
TDF  2.299  1.889  2.665  0.172  2.3  1.889  2.163  2.438 
TDFDDI  2.323  1.917  2.667  0.164  2.327  1.917  2.194  2.457 
ABC  2.459  1.733  2.698  0.126  2.485  1.733  2.404  2.544 
DDIABC  2.546  1.869  2.726  0.106  2.57  1.869  2.505  2.617 
DDI  2.746  2.675  2.78  0.02  2.747  2.675  2.732  2.762 
Figure 3 and Table 3 show that the probability distributions are almost uniform and ${\beta}_{i}\left({\varepsilon}_{X}^{i},{\varepsilon}_{Y}^{i}\right)$ values have considerable diversity and standard deviations among the viral strains. Hence, this observation means that even having point mutations can change the infection rates considerably and thus may lead to a need for more perfect adherence levels to the given therapy. Additionally, Fig. 3 implies that the initial viral strain of the patient plays a critical role in the progression of HIV dynamics. According to Table 3, NRTI therapy combinations yield 38.4% and 78% decrease in infection rate on average (among all therapies) (95% CI [36.2%–40.7%] and [69.7%–86.3%]) for the worst and best case scenario (having most and least resistant initial strain), respectively.
Table 3 ranks the possible NRTI combinations in terms of the resistance scores but ignores the side effects and costeffectiveness. Various sideeffects of NRTIs linked with mitochondrial toxicity (Holec et al., 2018). We present the possible sideeffects of the existing NRTIs in Table S3, and a detailed review can be found in the study of Montessori et al. (2004). The costeffectiveness of NRTI therapies is essential to maximize the expected survival times of the patients with minimized costs. Various mathematical models are available that compare treatments for costeffectiveness, and a detailed review of Mauskopf (2013) provides various essential results. Most of the models described in their study ignore the effect of drug resistance. Drug resistance is a crucial contributor to the expected costs. This study is only interested in the impact of drug resistance on the NRTI therapy outcomes, and we both ignore side effects and costeffectiveness.
Details of model simulations
In our simulations, we investigate the effect of the type of NRTI therapy, timing of the NRTI therapy, and adherence to the provided therapy on CD4+ T cell counts of the patients. All possible 21 mono and dual NRTI combinations of six inhibitors have been included in the simulations by considering their independent or additive actions. The initiation time of the NRTI therapy is considered within the first year after the patient became infected and denoted by τ. The adherence level of a patient to the provided therapy protocol is assigned to a real number α between 0 and 1, representing nonadherence to full adherence levels. After initiating the treatment with adherence level α in a day of the simulation, the patient takes drug(s) with probability α according to the parameters given in Table 1. Initial viral load, CD4 + T cell count, and macrophage cell count in the simulations are considered as 1 virion/ml, 10^{6} cell/ml and 150 cell/ml, respectively (HernandezVargas, 2019).
It is assumed that the patient is infected with one type of mutant strain with one to fivepoint mutations on the reverse transcriptase enzymes. In this way, five groups are constructed to include five different strains. These viral strains have been determined according to the frequency of presence in the Stanford HIV drug resistance database. These initial viral strains are denoted by G_{ij} where i = 1, 2, 3, 4, 5 denotes the number of the point mutations in the strain and j = 1, 2, 3, 4, 5 indexes the most frequently occurring examples in the dataset. We have performed our simulations with these 25 different initial viral strains having the following point mutations: ${G}_{11}=\left\{69D\right\},\phantom{\rule{0.25em}{0ex}}{G}_{12}=\left\{70E\right\},{G}_{13}=\left\{74I\right\}$, ${G}_{14}=\left\{151M\right\}$, ${G}_{15}=\left\{41L\right\}$, ${G}_{21}=\left\{69D,\phantom{\rule{0.25em}{0ex}}115F\right\},{G}_{22}=\left\{69D,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{23}=\left\{70R,215Y\right\}$, ${G}_{24}=\left\{67N,\phantom{\rule{0.25em}{0ex}}69D\right\},{G}_{25}=\left\{67N,\phantom{\rule{0.25em}{0ex}}70R\right\}$, ${G}_{31}=\left\{69D,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{32}=\left\{69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}115F\right\}$, ${G}_{33}=\left\{67N,\phantom{\rule{0.25em}{0ex}}69D,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{34}=\left\{67N,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{35}=\left\{67N,\phantom{\rule{0.25em}{0ex}}69D,\phantom{\rule{0.25em}{0ex}}70R\right\}$, ${G}_{41}=\left\{67N,69D,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{42}=\left\{67N,70R,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{43}=\left\{69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{44}=\left\{67N,69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}115F\right\}$, ${G}_{45}=\left\{65N,69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, G_{51} = {65N}, {69D, 70R, 115F, 215Y}, ${G}_{52}=\left\{69D,\phantom{\rule{0.25em}{0ex}}70R,74F,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{53}=\left\{41L,67N,69D,70R,\phantom{\rule{0.25em}{0ex}}215Y\right\}$, ${G}_{54}=\left\{65N,67N,69D,70R,215Y\right\}$, ${G}_{55}=\left\{67N,69D,\phantom{\rule{0.25em}{0ex}}70R,74I,\phantom{\rule{0.25em}{0ex}}215Y\right\}$. For instance, ${G}_{14}=\left\{151M\right\}$ strain has only one point mutation 151M and the rest of the amino acids are the same as wild type HIV1 virus.
Measuring the therapy success
It is essential to track the success of the given antiretroviral therapy by hindering the viral dynamics from the AIDS phase, i.e., by keeping the CD4 + T cell count as high as possible. The AIDS phase occurs when CD4 + T cell count is less than 200 cell/µl (Kitahata et al., 2009). Our primary criterion for the success of NRTI therapy is the occurrence and nonoccurrence of the AIDS phase after initiation of the therapy with some initiation timing τ and adherence level α, as was done in cohort studies (van Sighem et al., 2003). Note that it is also possible to increase the CD4 + T cell counts of patients during the AIDS phase by the initiation of the ARTs (Shoko & Chikobvu, 2019). However, here we are not analyzing what happens after the AIDS phase. Our primary goal is to determine how evolutionary dynamics under the NRTI therapy affect the occurrence of the AIDS phase.
All simulations start with one infected CD4 + T cell and one infected macrophage cell with one of the initial strains G_{ij}. The simulation final time t_{f} is considered 20 years, and therapy success/failure is determined according to the occurrence of the AIDS phase in 20 years. However, we note that the clinical goal of ART therapy is the full suppression of detectable viremia. In our simulations, total suppression of detectable viremia is equivalent to not developing AIDS after 20 years. However, the opposite is false: detectable (>200 copies/ml) suppression misses low copies of violent mutants, eventually leading to the AIDS phase. Therefore, we consider the AIDS occurrence as our output. In the clinic, therapy is redesigned if complete suppression is not observed. However, our simulations never redesign the treatment to distinguish between successful/failed drug combinations.
We run our simulations for randomly scattered 512 $\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\tau \right)\in \left[0,1\right]\times \left[0,365\right]$ pairs for predetermined initial strain G_{ij}. The success rate (SR) of a therapy is measured as the number of $\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\tau \right)$ pairs that lead to protection from the AIDS phase in all 512 $\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\tau \right)$ pairs. In Fig. 4, we show some representative simulation results of the multistrain withinhost model (5), starting with the ${G}_{51}=\left\{65N,\phantom{\rule{0.25em}{0ex}}69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$ strain under various mono and dual NRTI therapies with randomly scattered $\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\tau \right)$ pairs. For this simulation setup, nine out of 21 NRTI therapy protocols have considerable success in preventing the patient from the AIDS phase. The importance of adherence level (α) and initiation timing (τ) is evident from the figure for all cases. In some cases, such as the DDID4T combination shown in Fig. 4, the initiation timing considerably affects the success rates. Higher τ values yield therapy failure even at high adherence levels. As observed from the figure, the D4T3TC combination yields the best SR value by performing well for late initiation with perfect adherence levels. For the current case, the success of the D4T3TC combination is mainly due to the behaviour of the therapy in the higher initiation timing (τ) region.
While the importance of the adherence levels is evident from its direct relation with infection rates, the importance of the initiation timing is nonevident and should be explained here clearly. In Fig. 5, we illustrate the effect of initiation timing τ in our multistrain model (5) when initial strain and adherence level are selected as ${G}_{51}=\left\{65N,\phantom{\rule{0.25em}{0ex}}69D,\phantom{\rule{0.25em}{0ex}}70R,\phantom{\rule{0.25em}{0ex}}115F,\phantom{\rule{0.25em}{0ex}}215Y\right\}$ and α = 0.5. According to Figs. 5A–5B, τ = 50 yields successful therapy by maintaining healthy CD4 +T cell and macrophage cells at normal levels and declining the viral load to undetectable levels. On the other hand, when we assume the initiation timing as τ = 360, virologic failure and AIDS phase are observed in Figs. 5C–5D. According to our model (5), the main difference between early and late initiation timing is the diversity of viral strains at the initiation to therapy times. Late initiation to the therapy increases the probability of the occurrence of the more resistant strains, even if their ancestors are slowly growing. For example, as we compare Fig. 5B with Fig. 5D, the two generations of mutant strains occur when τ = 360 (Fig. 5D) while there exists only one generation of mutant strains when τ = 50 (Fig. 5B). The two generations of mutant strains yield viral rebound and failure of the therapy in Fig. 5D.
If we go back to Fig. 4, the NRTI combinations having boundary lines with relatively low slope values are more sensitive to increasing values of τ since these therapies yield high variance in IC_{50} values of possible viral strains mutated from the initial strain. Therefore, in our modelling framework, the late initiation is directly related to the variance of IC_{50} values corresponding to the initial strain and possible mutants. Thus, the level and type of the NRTI therapy should be planned so that the reoccurrence of the viral strains should be blocked depending on the initiation time τ. Additionally, in the reoccurrence phase of viral strains, nonperfect adherence to the therapy leads to the selection of resistant strains (Fig. 5D). In this case, two possible problems arise:

If the therapy protocol of the patient is updated, therapy is less likely to be successful than when therapy was first started.

The probability of infecting another person with more resistant strains increases, and the probability of having an AIDS phase increases for the infected person.
The existence of low viral loads of new mutated strains is enough for selecting these strains after antiretroviral therapy. Therefore, according to our simulations, initiation timing is as crucial as the adherence level to overcome the AIDS phase and to protect the possible susceptible persons from more dangerous scenarios.
The NRTI mutants are known to have epistasis effects, which implies that the viral fitness of the mutant strain depends on the existing genetic background. The epistasis effects may lead to the selection of diverse branches in mutant generations (Biswas et al., 2019). Epistasis of mutations can impact the values of IC_{50} and fitness costs. The data we used to train our IC_{50} values implicitly includes epistatic effects. The ANN model that predicts IC_{50} values for mutants is expected to learn the epistatic interactions. However, it is not completely unlikely that some unobserved data may have unpredictable epistasis. Nevertheless, that variant being underrepresented in the data implies its irrelevance in the clinic. On the other hand, the fitness costs of mutants are assumed to be fixed due to a lack of enough data. Nevertheless, as we explain later, this assumption should not significantly impact our claims.
Simulation results
Here we have simulated our multistrain withinhost model (5) for all possible initial strains G_{ij} to observe the effect of initial strains on success rates. All possible mono and dual NRTI therapies have been implemented for randomly scattered 512 $\left(\alpha ,\phantom{\rule{0.25em}{0ex}}\tau \right)\in \left[0,1\right]\times \left[0,365\right]$ pairs. The SR values of mono and dual NRTI therapies are calculated, and the wellperformed combination results are comparatively illustrated in Fig. 6.
In line with Fig. 3 and Table 3, the D4T3TC combination has been the best option for 20 out of 25 cases. The inhibitory potential of this combination is because of the pharmacokinetic parameters (see Table 1) of inhibitors, the drugresistance profiles of inhibitors (see Table 3), and their Blissindependent action on the target enzyme. Following the D4T3TC combination, the TDFD4T and D4TAZT combinations are observed to be in first place in four and one out of 25 cases, respectively. The strong relation between the infection rate of an initial strain (and possible new strains) and the corresponding success rate value is evident from the correlation between Figs. 3 and 6. For instance, according to Fig. 3, the D4T3TC combination yields fewer infection rates for most of the viral strains. Similarly, Fig. 6 shows that the D4T3TC combination has great success rates for most of the initial viral strains. We will later quantitatively analyze the relationship between the infection rates of the detected viral strains and the success rates of the given therapies.
According to our modelling framework, since the fitness cost of all strains is assumed to be the same, the initial strain is dominant when the patient is diagnosed. Moreover, as evident from Figs. 5B–5D, considerable mutational variations at low copy numbers exist besides the initial strain. However, only the dominant strain is likely to be detected (strains having higher than 200 copies/ml in blood Barletta, Edelman & Constantine, 2004) when a phenosense assay is implemented. Thus, the clinician would only observe the initial strain and maybe a few mutational variations (according to Figs. 5B–5D, only the initial strain can be observed when the patient is diagnosed) to decide on the NRTI therapy protocol. Therefore, it is inevitable to ask whether the only predictor of the success rate is the detected viral strains at the diagnosis.
The undetected viral strains play a vital role in estimating the success rate and finding an optimal therapy protocol—especially their infection rates. We have trained regression models that predict therapy outcomes based on the infection rates of the initial strain and its mutants—the mutants will be referred to as first, second, third, fourth, and fifth generations. The first generation is mutated from the initial strain, whereas the second is mutated from the first. For the regression model, we aimed to determine how many generations of the detected strain(s) should be considered to predict an optimal therapy. To answer this question quantitatively, we construct the ANN and MLR models for predicting the success rate of therapy from the infection rates of the existing mutant strains. We construct six ANN and MLR models denoted by G_{i} for i = 0, 1, …, 5.G_{i} denotes i − th generation of the detected strain(s) that has been considered in the inputs of the models. For instance, model G_{0} only assumes the infection rates of the detected viral strain(s), and model G_{3} considers the infection rates of the detected viral strain(s) and the first threegeneration mutants of this strain(s). In each generation of mutant strains, we use two values: mean and maximum values of the infection rates of the considered generation. Thus, together with the detected viral strain, the model G_{i} has 2i + 1 dimensional input. 2i input values denote the mean and maximum infection rates of i − th generation, and the remaining input value denotes the infection rate of the detected viral strain at the diagnosis. The graphical illustration of model G_{i} can be seen in Fig. 7.
Simulation results are given in Fig. 6 for 25 initial strains converted to the training data for the ANN and MLR models. 304 input–output relations have been obtained from various therapies having SR ≥ 0.02. For the ANN models, this data is divided into the train, test, and validation sets (70%, 15%, and 15%). Each G_{i} model having the ANN architecture is trained using the scaled conjugate gradient algorithm. Similarly, for the MLR models, 20% of the data is considered as a test set, and the remaining 80% is used in the training process. To test the prediction performances of the ANN and MLR models, we have generated an external test dataset by simulating the model (5) with 25 random initial strains having onetofivepoint mutations, and 314 test sample is obtained. Additionally, to observe how well our ANN and MLR models classify the therapies as successful (SR > = 0.5) and unsuccessful (SR < 0.5), the area under the receiving operating curves is measured for both the ANN and MLR models.
We illustrate the regression and classification performances of the ANN models on the training and test sets in Fig. 8. Figure 9 shows similar predictive performance metrics of the MLR models on both the training and test sets. The mean square error (MSE), linear correlation coefficient (R), and area under the curve (AUC) metrics are presented for six G_{i} models having the ANN and MLR architectures. According to the test set performance of the models, model G_{2} gives better MSE, R, and AUC values with both the ANN and MLR architectures. That means considering the infection rates of both the detected strains and the first two mutant generations of the detected strains led to better predictions.
On the other hand, the G_{0} type models yield relatively poor regression and classification performances, i.e., considering only the infection rate of the detected strains is not enough to estimate better therapy protocols. This implies that the possible undetected mutant generations should also be taken into account in determining the therapy protocols. Nevertheless, there is a threshold on the number of mutant generations that must be considered. Figure 8 (for ANN architectures) and Fig. 9 (for MLR architectures) show that models G_{3}, G_{4} and G_{5} overfit the data and yield less accurate predictions than the model G_{2} for both architectures. Additionally, for each G_{i} model, the ANN architecture yields a better approximation for the SR values than the MLR architectures.
Discussions and Conclusions
In this study, we have proposed a multistrain withinhost model of HIV infection with timedependent NRTI therapy. Drugresistant strains have been assumed to initiate the infection for the patients, and six available NRTI inhibitors with mono and dual combinations have been implemented in the simulations for various initiation timing and adherence levels. To assess the drug response curves with the IC_{50} values of the NRTIresistant strains, artificial neural network models are trained for each inhibitor by using the Stanford HIV drug resistance database. To describe timedrug efficiency and timeinfection rate curves, pharmacokinetic parameters of the inhibitors have been calculated and hybridized with the corresponding IC_{50} values. We have designed our simulation environment to determine the effect of initial strains, initiation timing for the therapy protocol, and adherence levels to the given drug usage schedule on the occurrence of the AIDS phase within 20 years after infection.
According to our modelling framework, the success rate of the NRTI therapies in case of late initiation has led to the availability of more resistant viral strains, and then the resistant strains become dominant in the host plasma after an initial decline of the detected strain. Although some mathematical models assume implicitly that the initiation timing does not affect the successfailure of the therapy (Dixit & Perelson, 2004; Rong, Feng & Perelson, 2007), our multistrain model catches the penalty of late initiation since the late initiation was proven to block the therapy success in various experimental results (Kitahata et al., 2009; van Sighem et al., 2003). Our simulation results have shown that in the case of the late initiation to therapy, the efficiency of the therapy should be far more than the early initiation case to prevent the possible AIDS phase.
We have shown that D4T3TC, D4TAZT, and TDFD4T combinations are less likely to result in treatment failure. These inhibitors have been seen to provide fewer infection rates due to their pharmacokinetic parameters and IC_{50} values in the presence of various viral strains. According to our results, the success rate of accurately predicting the best therapy depended on the composition of detected strains and their possible further mutants. This observation implies that the emergence of new mutants from the initial strain is likely to have a considerable effect on the success of the therapy. Thus, it is more reasonable to suggest the optimal therapy combinations for the patients by considering the detected viral strain and the undetected mutant, which most likely were generated from the detected strain.
The most important message of this article is that the undetected viral strains, at the diagnosis, may have considerable effects on therapy outcomes. Specifically, double mutants of the detected viral strain should be taken into account even if they were not detected. Earlier studies, such as Stanford HIVdb (Talbot et al., 2010), HIVgrade (Obermeier et al., 2012), REGA (Van Laethem et al., 2002), and ANRS (Agence Nationale de Recherches sur le SIDA, Meynard et al. (2002)) predicted the best possible therapy protocol. REGA is a rulebased model and was developed by scientists at Rega Institute for Medical Research and University Hospitals, and classifies the isolates as susceptible, intermediate, and resistant (Van Laethem et al., 2002). ANRS Meynard et al. (2002); Singh (2017) is also a rulebased computational resistance classifier based on a linear combination of mutations. However, the undetected viral strains may lower the prediction power of such models. We have shown that a multistrain withinhost model (5) can help estimate undetected mutant strains and their role in optimal therapy selection.
A possible criticism of our model is that each mutant strain should have a unique fitness cost. However, we assume a constant factor for all mutants. To our best knowledge, there is not much data for specific strains to construct a machinelearning model as we did for the IC_{50} values. According to the theory, fitness costs can play a role in selecting resistant strains, which can alter our success rate. However, the fitness costs would affect the dynamics more at low drug concentrations. Luckily, the phase changes (AIDS or no AIDS) occur at relatively high adherence levels, which implies a relatively high concentration.
Our modeled treatments include only NRTIs, but current clinical practice includes additional drugs (Aguilar et al., 2022). Indeed, including the other components of ART would add to the realism. However, it is known that different classes of HIV drugs generally interact independently (Rosenbloom et al., 2012; Jilek et al., 2012). By the independence assumption, the relative ranking of NRTI therapies is relevant to consideration for ART. However, we would like to openly indicate that our model is not designed to suggest a better first line of treatment but rather to relatively rank NRTI combinations in a multiscale model.
This study has investigated the effect of NRTI inhibitors, which are the most important members of Highly Active Antiretroviral Therapy (HAART) (Achhra & Boyd, 2013). Since the Stanford drug resistance database also includes the genotypephenotype data of protease inhibitors (PI), nonnucleotide reverse transcriptase inhibitors (NNRTI), and integrase inhibitors (II), some future studies may include these groups of inhibitors with possible mono, dual or triple drug combinations. Some existing HAART protocols may also be simulated through such a modelling framework. On the other hand, we have not considered the toolate initiation of the NRTI therapy at considerably low CD4 + T cell levels because of the failure of simulated therapy protocols in such situations. Some future works may also investigate more comprehensive therapies to prevent patients from the AIDS phase when they are diagnosed too late.
Supplemental Information
The ANN hyperparameter tuning results
Mean square error (MSE) and linear correlation coefficient (R) for artificial neural network (ANN) models with different topologies (number of nodes per layer) for predicting logarithmic foldchange of IC50 values for the six NRTIs. Genotypephenotype data of each NRTI is split as external test set and internal trainvalidationtest sets with 15% and 85%. Internal data sets are divided into training, validation, and test sets (80%, 10%, and %10). Internal data set is used for training and model selection. In the model selection process, we trained 20 models four times, and every 20 models have been tested on the internal test sets, and the models giving the best MSE scores are selected for each NRTI. Thus, four models are averaged to get the final model for each inhibitor. The selected models are tested on the external test set, and the resulting performance metrics are presented below. According to the table, two topologies yield relatively better results: the model has one layer with five neurons, and the model has two layers and three neurons in each layer. To make the model under consideration simple enough, we chose the model with a fiveneuron hidden layer.
Comparison of the ANN and MLR models
Mean square error (MSE), linear correlation coefficient (R) and area under the curve (AUC) metric values for linear regression (LR) and artificial neural network (ANN) models for predicting foldchange values in IC50 values of the six NRTIs.
Pharmacokinetic parameters of the nucleotide reverse transcriptase inhibitors with considerable side effects
ROCs and corresponding AUC values for the ANN models
The ANN models are used to classify the given strains as resistant (Fold Change ≥ 3) and susceptible (Fold Change<3). The corresponding receiver operating characteristic curves with area under the curve (AUC) values are presented for each NRTI.
Regression performance of the six MLR models for each NRTI to predict logarithmic fold change values (log(FC)) of the mutant strains existing in the data
The xaxis of the figures denotes logarithmic fold change value for all existing mutant strains in the data and yaxis denotes corresponding predictions of the MLR models. For each MLR model, linear correlation coefficient (R) and mean square error (MSE) metrics are specified to measure the ability of these models to fit the existing real data.
ROCs and corresponding AUC values for the MLR models
The MLR models are used to classify the given strains as resistant (Fold Change ≥ 3) and susceptible (Fold Change<3). The corresponding receiver operating characteristic curves with area under the curve (AUC) values are presented for each NRTI.