This proposal was reviewed and approved by the University of Rochester Human Subjects Review Board (protocol number RSRB00068884). Consent was waived by the review board due to de-identification of the data set. The analysis in this paper is presented in compliance with Center for Medicare Services (CMS) current cell size suppression policy ²⁰ . Data were coded such that patients could not be identified directly in compliance with the Department of Health and Human Services Regulations for the Protection of Human Subjects (45 CFR 46.101(b)(4)).

We obtained medical encounter data from all patients with an HIV diagnosis in the University of Rochester Medical Center’s electronic medical record system (EMR) between 2011–2016, including age, gender, race, ethnicity, and VL measurements. There were a total of 1,892 patients with at least one VL measured, with 1,576 of these patients having at least three VL measurements.

Measurements ≤48 copies/mL, present as categorical values “NEG", “POS < 20", or “POS < 48" were transformed into numerical values of 0, 20, and 48 respectively. The deidentified study data containing only viral load measurements and relative time are available at https://doi.org/10.5281/zenodo.1313245 ²¹ .

Analyses were performed on a Windows 8 server with Intel(R) Xeon(R) CPUs E5-2620 v2 @ 2.10GHz and 256GB of RAM. Python 2.7 was used for most data mining and machine learning under Spyder v.3 installed from Anaconda2 (64-bit). The default packages available in Anaconda were used for analysis, including, but not limited to: NumPy, scikit-learn, SciPy, datetime, csv, math, Matplotlib, pip, operator, copy, random, and time. Using pip we installed the webcolors and pydotplus packages for rendering a decision tree. SQLite was used to store, query, and clean ~the data. Analytic code is available for download at https://github.com/SamirRCHI/Viral_Load_Data_Categorization.

Since VL data is asynchronous and noisy, with variable numbers of data points for each subject, we excluded patients with ≤ 2 VL measurements as too few to accurately assess VL patterns. Based on temporal patterns of VL described in the literature, the VL pair distribution of our data ( Figure S1), and a further extensive investigation into the data, we hypothesized six potential temporal VL patterns, defined in Table 1 and illustrated in Figure 1.

It is important to note that these definitions are pattern based, and do not explicitly select absolute VL cutoff levels or a specific temporal window, as other reports have done ^{8– 10, 13} . This has the advantage of allowing the absolute VL levels and critical time windows to emerge from the analysis. It also does not preclude incorporation of absolute levels (e.g. VLM >400) at a later stage into the pattern specification.

Mathematical notations for this work are described in Table 2. We next designed a feature vector to capture characteristics that would allow us to distinguish between VL patterns. VL values at the lower limit of detection are a function of the specific assay used, and appear in our data set as 0, 20, and 48 copies/mL ( Figure S1). Thus, plots of the log ₁₀ transformed data have discretely spaced values at the lower level of detection, capturing the undetectable range of viral load. Additionally, we adjusted the data by log ₁₀[ V L + 10] to avoid log ₁₀[0]. The addition of 10 to VL (instead of 1) is used to minimize the distance between the undetectable values: 0, 20, and 48 (copies/mL). Thus, in our notation, all the values related to viral load are assumed to have been adjusted to this measure. For example, min _{V L} = log ₁₀[0 + 10] = 1 and max _{V L} = log ₁₀[10 ⁷ + 10] ≈ 7.

Using the transformed VL data, we next extract several relevant features of the VL measurements over time. These features are used for machine learning classification of individual patient VL time series, and designed to distinguish patterns in VL change while minimizing the effects of noise. We do not limit feature extraction based on the total elapsed time of viral load measurements because the optimal time-point for determining viral load class is not well established. The attributes for feature extraction are: relative area of viral exposure, weighted recency reliability, adjusted maximal difference, and interquartile range. The definitions include:

A ˙ p = ∑ i = 2 V L M p V L → p , i + V L → p , i − 1 2 ( t → p , i − t → p , i − 1 ) ( m a x V L − m i n V L ) ( t → p , V L M p − t → p , 1 ) ( 1 )

w e i g h t → p , i = 1 t → p , V L M p − t → p , i + 1 ( 2 )

w R p = w e i g h t → p • V L → p ∑ i = 1 V L M p w e i g h t → p , i ( 3 )

d e v → p , i = | w R p − V L → p , i | ( 4 )

w R R p = 1 median( d e v → p ) + 1 ( 5 )

grnd ( x ) = { − 1 0 x < 0 x ≥ 0 ( 6 )

D ˅ p = grnd ( w R p − V L → p , argmax ⁡ ( d e v → p ) ) ⋅ max ⁡ ( d e v → p ) ⋅ w R R p ( 7 )

I Q R p = Q 3 ( V L → p ) − Q 1 ( V L → p ) ( 8 )

Machine learning methods for cluster classification were compared by calculating F ₁ scores, the harmonic mean of precision and recall ²² , defined by Equation 9– Equation 11.

p r e c i s i o n = T r u e P o s i t i v e T r u e P o s i t i v e + F a l s e P o s i t i v e ( 9 )

r e c a l l = T r u e P o s i t i v e P o s i t i v e ( 10 )

F 1 = 2 ∗ p r e c i s i o n ∗ r e c a l l p r e c i s i o n + r e c a l l ( 11 )

Here we formally define keywords appearing in the analysis: Let Feature extraction be the process of determining the values Ȧ, wRR, D ˅ , and IQR from a set of patients (using their viral load patterns) with the formulations given above. Then a feature vector ( F → _p ) contains the values Ȧ _p , wRR _p , D ˅ p , and IQR _p extracted from patient p’s viral load pattern. The words sample or point are also used here RVL (black; n= 237) and HVLS (purple; n=316) clusters. interchangeably. The term feature ( F) can be thought of as a column vector for all patients in the dataset consisting of the four attributes: F _Ȧ , F _wRR , F D ˅ , and F _IQR . Finally, the terms label assignment, VL pattern membership assignment, patient categorization, and prediction, all refer to the same principle: To assign the most appropriate label which characterizes the viral load pattern of a patient. However, while the principle is the same, the method of assigning such an appropriate label differs depending on the categorization or the learning method used.

We began by transforming viral load data by min-max normalization ²² to equally weight the temporal features of the VL series ( Equation 12). That is, we normalize the features, F, to a range between [0, 1] using Equation 12 where F ^* = f( F).

F ∗ = f ( F ) = F − min ⁡ F max ⁡ F − min ⁡ F ( 12 )

Next, we examined each of the four features for all patients with ≥ 3 viral load measurements ( N = 1,576 patients), and did not find distinct bi-variate clustering ( Figure S2). A feature correlation coefficient analysis ( Supplementary Table S1) revealed that the Adj MD feature is linearly independent of Area and wRR. In contrast, there is modest linear dependence between IQR and Adj MD, and between Area and both wRR and IQR. As expected, the largest linear dependency is between wRR and IQR. These results suggest the separation between viral load patterns will be most noticeable between the Area and the Adj MD features - as we designed them to be. Also, although Adj MD is dependent upon wRR, we find that their correlation coefficient is very low (0.033).

We then performed hierarchical clustering of the individual subject VL patterns using a Euclidean distance metric and Ward’s linkage criterion ²³ to minimize the total within-cluster variance. Patients showed a clear separation into 5 distinct groups, which had clinical significance ( Figure 2 and Figure 3). The cluster with the lowest viral loads and the highest weighted recency reliability (n=442) corresponds to the DSVL patient group. The patients corresponding to the SHVL group (orange; n=46) exhibited the highest relative Area and very low IQR. Compared to the DSVL cluster, the blue cluster (n=535) has slightly greater area and IQR with a significant difference in the weighted recency reliability. Using this information, along with the general patterns shown by Figure 3, we identify this as the SLVL group. The algorithm also identifies the RVL (black; n=237) and HVLS (purple; n=316) clusters. The RVL cluster has a low weighted recency reliability and high IQR. In contrast, the HVLS cluster has a lower area, higher weighted recency reliability, indicating little variation in the terminal portion of the VL time series, and most importantly very low adjusted maximal differences ( Figure 4).

VL patterns are similar within clusters, and dissimilar between clusters Figure 3. Interestingly, there are patients within each cluster whose last VL measurement occurs near 1,827 days. This is equivalent to the full span of five years of VL monitoring data set. This suggests that these clusters don’t disappear after some elapsed time, but rather each type of pattern can be found at virtually any time point.

We found large VL spikes within the time series of the HVLS group. We hypothesize that this may be due to the asynchronous timing of measurements between subjects, the natural variation in biological responses, or patient variability in adherence to therapy. This observation also reflects one limitation of asynchronous outcomes data sampling, which lacks a “completion" endpoint characteristic of most prospective, randomized clinical trials. If measurements ended at a spike, the adjusted maximal difference feature may be weighted in the favor of the patient being classified as RVL. This may indicate that some patients classified as suppressing their viral loads should have been classified as having rebounding viral loads. Alternatively, may indicate that these features do not restrict a patient to forever to one category, but allow for dynamic classification as a function of biological or therapeutic responses.

Using the same data set, we next compared our VL pattern categorization method to those previously published in the literature. Visually, we find that the SLVL group detected by our method is very similar to the LLVR group defined by Greub et al. ( Figure 5). Furthermore, it appears that the methods trying to capture SHVL, viral rebound, and viral failure patients did not succeed as well as the identification of SHVL and RVL patients in our method. RMVL repeat continuous visually appears to have performed very well in identifying patients whom have suppressed their viral loads. However, the results suggest that our analysis performs slightly better in identifying the suppression group (HVLS), as we find that the last VL measurements (black dots in Figure 5) are consistently low using our method.

The other methods may not have performed as well as they rely on a window or a consecutive pair measure, which may be too subjective for assigning VL pattern membership. Furthermore, notice that patients with baseline VL<200 ( Figure 5) contain VL patterns which can reach as high as 10 ⁶ copies/ml, which is in contrast to Rose et al.’s assumption that these patients have consistently low viral variation. Lastly we wish to emphasize that while some of these categorization methods are successful in identifying a specific group of patients, our method is unique as it attempts to associate each VL pattern to a specific category, without using categories such as “Not Suppressed”, “Unspecified”, or “Omitted”.

We next used the classes identified by hierarchical clustering to compare several machine learning models, with the goal of identifying methods that could be trained to prospectively assign HIV patients to VL categories (i.e. SHVL, SVL, SLVL, DSVL, and RVL). Unsupervised learning (e.g. hierarchical clustering) is useful for establishing the data structure of VL categories and their locations in the feature vector space. Once the model is established (e.g. cluster boundaries), supervised learning methods are better suited for prospective cluster assignment, given a robust "ground truth" for model training, as they do not depend on re-analysis of the entire population.

To this end, we compared the predictive power of several supervised learning methods for HIV cluster assignment, including: k-nearest neighbors (kNN), decision tree, support vector machine (SVM), Adaboost, and random forests. Models were trained on the original data set, and we then ranked their prediction power by their average F ₁ score derived by leave-one-out cross-validation (LOOCV) on the clustered results ( Table 3). We compared the ability of these methods to reconstruct the originally identified clusters, even when allowing for variability in cluster numbers (e.g. kNN with k={7, 9}, or DT without a maximum depth specification). All methods performed comparably well, with the notable exception of Adaboost. This generally high performance was expected because the VL pattern categories are well-separated as a result of the clustering. k-Nearest Neighbors and k=5, was computationally efficient and yielded the best results in Table 3.

We next considered the trade-off of predictive precision versus model interpretability. Critical clinical evaluation of machine learning results is important to protect against mis-categorization and clinical error. For this reason, many have advocated using models that are more clinically interpretable. kNN is dependent upon the entire training set for prediction, as it does not inherently “learn” patterns ²⁴ , hence it does not meet our interpretability criteria. In comparison, SVM offers a simpler model, but it’s results could be non-intuitive for clinicians. And although Decision Trees offer the best interpretability, overly complex trees may be generated, as occurred in our study ( Figure S3).

We found that pruned decision tree rules, with a maximum depth of 5 levels, met this interpretable criteria, however at a slight cost to the predictive power ( Table 3). The extracted decision rules are shown in Figure 6. Each category has a rule with a high proportion of true positive samples following the rule relative to all samples for the category (support). Similarly, a high proportion of the predicted class was found in the rule (precision), indicating that the rules can be summarized into a majority rule. Note that the sum of the support does not necessarily add up to one for each class because some samples belonging to that class may have been otherwise placed into a different rule, making the precision of that rule weaker.

As an alternative interpretable model we explored the use of centroid cluster summarization, which is often used in clustering algorithms, and is flexible enough to accommodate different centroid determination methods ^{22, 25} . To determine the effects of different centroid determination algorithms, we compared seven different methods: multidimensional mean, multidimensional median, best representative center, bounding box method, smallest disk method, polyhedral center, and a novel “push and pull" (PnP) method inspired by force-directed graph drawing such as the Fruchterman-Reingold’s algorithm ^{26, 27} (see Supplementary File 1). Force directed clustering methods maximize inter-cluster center distances, while minimizing intra-cluster distance, and are the basis for modularity clustering in graph theory ²⁸ .

We then combined the centroid cluster summarization approach with a radius-based classification prediction algorithm. Let c _i be the ith cluster center with corresponding radius r _i , where r _i is calculated as the distance to the farthest intra-cluster sample from c _i , then for a new sample s choose its predicted cluster membership j such that ‖ s − c j ‖ 2 r j is a minimum. We refer to this method as radial normalization classification.

Comparing the representative F ₁ power of the centroid radial normalization methods (italicized in Table 3) to common machine learning algorithms, we find that the centroid interpretation loses some predictive power. However, the centroid summary is highly interpretable because the entire model can be expressed concisely ( Supplementary Table S2), and understood clearly. For example, a clinician classifying a patient by VL time series values would compare observed feature values with the ranges given in Figure 6, and find which classification the patient’s data fits best within. In the case of the centroid method, if an observed value appears to fall in multiple categories, then they should be assigned to the one closest to the center (this allows a clinician to cross-check model predictions).

HIV patient viral load states are often fluid, with class changes (e.g. SHVL → HVLS) occurring due to therapy, viral genetics, social and other factors. To examine this aspect of classes, we use the k-Nearest Neighbors (k=5) model, fit to the original clusters, to predict the class state of each patient with ≥ 3 VL measurements using only partially retained VL data. For example if a patient has 6 viral load measurements, then we predict the class state at 3, 4, 5, and 6 VLM, which may yield SHVL → SHVL → RVL → HVLS as its prediction. We then constructed a state-transfer network using the trace-route method ²⁹ , revealing several interesting relationships: