Whole‐brain computational modeling reveals disruption of microscale brain dynamics in HIV infected individuals

Abstract MRI‐based neuroimaging techniques have been used to investigate brain injury associated with HIV‐infection. Whole‐brain cortical mean‐field dynamic modeling provides a way to integrate structural and functional imaging outcomes, allowing investigation of microscale brain dynamics. In this study, we adopted the relaxed mean‐field dynamic modeling to investigate structural and functional connectivity in 42 HIV‐infected subjects before and after 12‐week of combination antiretroviral therapy (cART) and compared them with 46 age‐matched healthy subjects. Microscale brain dynamics were modeled by a set of parameters including two region‐specific microscale brain properties, recurrent connection strengths, and subcortical inputs. We also analyzed the relationship between the model parameters (i.e., the recurrent connection and subcortical inputs) and functional network topological characterizations, including smallworldness, clustering coefficient, and network efficiency. The results show that untreated HIV‐infected individuals have disrupted local brain dynamics that in part correlate with network topological measurements. Notably, after 12 weeks of cART, both the microscale brain dynamics and the network topological measurements improved and were closer to those in the healthy brain. This was also associated with improved cognitive performance, suggesting that improvement in local brain dynamics translates into clinical improvement.


| INTRODUCTION
There were approximately 37.9 million people globally living with the HIV in 2018 according to World Health Organization (WHO, https:// www.who.int/gho/hiv/en/). Since the introduction of antiretroviral therapy (ART), there has been a significant decrease in the mortality rate of people infected with HIV, and a dramatic decrease in the incidence of HIV-associated dementia (Saylor et al., 2016). However, the prevalence of HIV-associated neurocognitive disorders (HAND) is increasing (Nabha, Duong, & Timpone, 2013). The cause of the high prevalence of HAND remains unclear. It is likely multifactorial, including early injury prior to starting antiretroviral drugs (HIV infects the brain shortly after transmission), chronic mild neuroinflammation and possibly neurotoxicity of antiretroviral drugs (Chang & Shukla, 2018;Fois & Brew, 2015).
Overall, there is evidence that brain injury is present and measurable via MRI in HIV infected individuals. However, most of the previous studies have investigated one modality at a time, whereas the integration of multimodalities in HIV-related studies has not been fully explored yet. Large-scale whole-brain dynamic modeling is a promising approach to further quantify the integrated contribution of structural and FC in CNS injury. This approach allows simulating resting-state fluctuations emerging from the interaction between brain regions, constrained by the anatomical connections derived from DTI, and effectively integrates FC and structural connectivity (SC) (Honey et al., 2009;Mollink et al., 2019;Surampudi et al., 2019;Wang et al., 2019). Furthermore, this approach has already provided some insights in connectome disruption in other neurological disorders such as Alzheimer's disease (AD) and Parkinson's disease (Alderson et al., 2018;Deco & Kringelbach, 2014;Demirtas et al., 2017;Jirsa, Sporns, Breakspear, Deco, & McIntosh, 2010;Proix, Bartolomei, Guye, & Jirsa, 2017). A novel whole-brain modeling technique, named relaxed mean field dynamic modeling (rMFM) (Wang et al., 2019) has been proposed to simulate local brain dynamics. Previous wholebrain modeling studies assumed that local microscale properties, the recurrent connection strengths and subcortical inputs, were the same across entire brain, while the rMFM modeling method relaxed these two parameters to be heterogeneous across different brain regions.
Two microscale brain properties, recurrent connection strengths and subcortical inputs, can be derived from this generative brain dynamic modeling by tuning the model to fit the simulated FC to empirical FC.
In this study, we applied the rMFM (Wang et al., 2019) to assess CNS changes in a cohort of HIV infected treatment-naïve patients at baseline (HIV + BSL) and after 12 weeks of combination antiretroviral therapy (cART) treatment (HIV + 12 wk), and compared them to HCs.
The overreaching goal was to investigate potentially new imaging biomarkers that could be sensitive to CNS changes and thus helpful in monitoring CNS disease progression and response to treatment. First, we built the rMFM whole-brain dynamic models for the three groups (HC, HIV + BSL, HIV + 12 wk), respectively. We then compared the microscale brain properties (the recurrent connection strength and subcortical inputs) across the three groups. This was followed by investigating the topological changes among the three groups using traditional graph theoretical analysis from both global and regional perspectives. Subsequently, we investigated the association of nodal graph theoretical measurements with microscale brain properties and neuropsychological tests scores.

| MRI data collection
MRI was performed on a 3 T Siemens MAGNETOM Trio MRI scanner (Siemens Healthineer, Germany) equipped with a 32-channel head coil. T1-weighted three-dimensional magnetization-prepared rapid acquisition gradient echo images were acquired, with repetition time (TR)/inversion time (TI)/echo time (TE) = 2,530/1,100/3.44 ms, voxel size = 1.0 × 1.0 × 1.0 mm 3 , flip angle = 7 , bandwidth = 190 Hz/pixel. Diffusion weighted imaging (DWI) data were acquired with the following parameters: TR/TE = 34.85/7.12 ms; 10 b = 0 s/mm 2 images; 60 diffusion weighting images with b = 1,000 s/mm 2 and direction uniformly distributed on a unit sphere; voxel size = 2 × 2 × 2 mm 3 , bandwidth = 1,502 Hz/pixel. A double-echo gradient echo field map sequence was acquired with the same resolution as the DWI sequence and was used to correct for distortion caused by B0 inhomogeneity. Resting-state fMRI data were acquired using a gradient echo-planar imaging sequence, with TR/TE = 2,000/30 ms, voxel size = 4 × 4 × 4 mm 3 , 150 time points, flip angle = 90 , bandwidth = 1,562 Hz/pixel. During the entire 5-min resting-state fMRI series, participants were instructed to keep their eyes closed and avoid falling asleep.

| Overview of data processing
The processing pipeline schematically shown in Figure 1 has the following steps: 1. We preprocessed the T1-weighted (T1w) image, DWI, and rsfMRI, and constructed the empirical SC and FC (Section 2.3).
2. The rMFM whole-brain model was built for each group (Section 2.4): We split the empirical SC and FC into training and test datasets for each cohort. We then derived the rMFM model parameters using the training dataset and validated each rMFM model using the test datasets (Section 2.5). The microscale brain properties, namely recurrent connection strength, denoted as w, and subcortical input strength, denoted as I, were derived from each group were then compared among HIV + BSL, HIV + 12 wk and HC subjects.
3. Graph topological analysis on empirical FC and SC (Section 2.6): We investigated the topological changes in empirical FC and correlated the recurrent connection strength w, and subcortical inputs I with the nodal graph theoretical measurements.
4. The rMFM modeling and graph theoretical analysis were validated in three steps (see supplementary material): (a) split the data into different training and testing groups, (b) doubled the rMFM optimization steps, and (c) reproduced our results using a finer segmented brain atlas, Destrieux atlas (Destrieux, Fischl, Dale, & Halgren, 2010), which includes 148 cortical regions. 5. Neuropsychological assessment scores were compared among cohorts, and their relations with graph theoretical measurements were also investigated (Section 2.7).

| Anatomical T1w image analysis
Three-dimensional high resolution T1w images were brain extracted using Brain Extraction Tool (Smith, 2002) to remove nonbrain tissue.

| DTI data processing and whole-brain tractography
Preprocessing of DTI data was reported in detail previously (Zhuang et al., 2017). Briefly, the b0 images and diffusion-weighted images were F I G U R E 1 Overview of processing pipeline. (a) We collected diffusion weighted image (DWI), T1-weighted (T1w) structural image, and resting-state functional MRI (rsfMRI) for each HIV+ and healthy control (HC) subject. (b) For each subject, the T1w image was segmented using Freesurfer. DWI and rsfMRI were preprocessed and coregistered to T1w image. Whole-brain tractography for each subject was derived from DWI image, and a structural connectivity (SC) matrix was generated for each subject. The streamline count for each pair of brain regions constitutes the SC matrix. Regional BOLD time series were extracted for each subject. The Pearson's correlation coefficient of the BOLD time series of two brain regions constitutes the functional connectivity (FC) matrix. (c) Whole-brain dynamic modeling, specifically, relaxed mean field dynamic modeling (rMFM), is used to study the neuronal dynamic changes between HIV+ subjects and HC subjects. (d) We also used graph theoretical analysis to investigate the topological changes by comparing the SC and FC for HIV+ subjects with HC subjects motion corrected using a 6-degrees of freedom rigid-body registration and field maps were used to correct the susceptibility induced distortion, using FUGUE in FSL (Jenkinson, 2003). DTI processing and structural connectome construction were then performed using the population-based structural connectome pipeline (Zhang et al., 2018). A reproducible probabilistic whole-brain tractography algorithm (Girard, Whittingstall, Deriche, & Descoteaux, 2014;Maier-Hein et al., 2017) was used to reconstruct streamlines. The tissue partial volume estimation maps obtained from the anatomical T1w image helped reduce the tractography bias. The particle filtering tractography was used to reconstruct more reliable streamlines (Girard et al., 2014). After whole-brain tractography was generated, the Desikan-Killiany atlas and Destrieux atlas were used to construct the structural connectome represented by 68 × 68 matrix and 148 × 148 matrix, respectively. The individual SC matrix was used for graphic theoretical analysis. Subjects were randomly split into training and test groups, controlling for age, so that there was no significant age difference between the training and testing group, and between the three cohorts (see Supplemental Table S1). For the whole-brain dynamic modeling, the SC matrices were averaged across subjects in training and test datasets, respectively.
Detailed preprocessed steps are included in the supplementary materials. The preprocessed fMRI data were then further cleaned by using FIX Salimi-Khorshidi et al., 2014), an ICA-based noise removal software that automatically classifies noise components of the resting-state functional data. Head motion, white matter, and cerebrospinal fluid time series were regressed out as nuisance regressors. Two HIV + BSL subjects and one HC subject were removed from rsfMRI analysis due to head motion greater than 2.5 mm. We used a two-sample t test to ensure there was no significant head motion difference between HIV+ and HC cohort (HIV+ mean displacement = 0.748 mm, HC mean displacement = 0.645 mm, p = .065).
Mean BOLD signals were then extracted from each parcel in the Desikan-Killiany and Destrieux atlases. Pearson's correlation of the mean BOLD signals between each pair of parcels was calculated yielding an FC matrix for each subject. The FC matrix was then Fisher's r-to-z transformed, resulting in a zFC matrix. The zFC matrices for all subjects were used for graph theoretical analysis. In preparation for rMFM whole-brain dynamic modeling, the zFC matrices were averaged across subjects within training or test dataset, for each age-controlled cohort (HIV + BSL, HIV + 12 wk, and HC) respectively (Supplemental Table S1).

| Simulated FC using rMFM model
The whole-brain dynamic mean-field model (MFM) (Deco et al., 2013; has been widely used to derive spontaneous brain activity from SC (Deco et al., 2018). Here, we used a modified version of MFM, the rMFM (Wang et al., 2019), which assumes the recurrent connection strengths w and subcortical inputs I are not uniformly distributed in the brain. This model has been previously proved to improve the FC simulation by 53% over the original MFM (Wang et al., 2019). Using rMFM, we simulated the neural activity for each cortical region, and then used the Balloon-Windkessel hemodynamic model (Buxton, Wong, & Frank, 1998;Friston, Harrison, & Penny, 2003;Friston, Mechelli, Turner, & Price, 2000) to convert neural activity to simulated BOLD signal. The simulated FC was then calculated using Pearson's correlation of the BOLD signal for pairwise cortical regions. More details of the rMFM model and its mathematical relations to various parameters are given in supplementary materials.

| Model parameters estimation
Empirical SC and FC for each group (HIV + BSL, HIV + 12 wk, and HC) in the training dataset was used to estimate the model parameters. There are 138 parameters to be optimized when we use the Desikan atlas (68 recurrent connection strength w i , 68 subcortical input strength I i , global scaling factor G, and noise coefficient σ).
The optimization steps were based on the expectationmaximization algorithm in dynamic causal modeling (Friston et al., 2003), and were detailed in (Wang et al., 2019). In this study, 500 iterations were performed for each run, and the final optimum estimated parameters for the rMFM model were chosen from the highest Pearson's correlation coefficient between the empirical FC and simulated FC in the training dataset from the 500 iterations (highest similarity for each cohort in Supplemental Figure S1a). The corresponding 138 model parameters were stored, shown as one column in Supplemental Figure S2.
We noticed that the random initialization parameters had a small effect on the final model parameters, so we repeated the entire optimization process with 25 different random initializations for each cohort. We then chose the top five sets of model parameters, that is, the model parameters corresponding to the five highest Pearson's correlation coefficients between empirical FC and simulated FC (the five highest similarities in Supplemental Figure S1b).
We then fed the empirical SC in the test dataset for each group to its own rMFM model, to validate these three rMFM models by calculating the similarity between the simulated FC and empirical FC for each group. After each model for HIV + BSL, for HIV + 12 wk, and for HC was validated, we obtained three rMFM generative models representing the three cohorts. The recurrent connection w = {w 1 , …, w n } and subcortical inputs I = {I 1 , …, I n }, where n = 68 for Desikan-Killiany atlas, in each rMFM model represent the microscale brain dynamics for each group and were compared among groups.
To ensure 500 iterations for each optimization run were sufficient, we also repeated the process with 1,000 iterations, with 15 different random initializations for each group. In addition, we repeated the rMFM modeling process using a different atlas to validate our results. To do this, we constructed the empirical SC and FC using Destrieux atlas, which yields 298 model parameters to be optimized.
Detailed results from these analyses are reported in the supplementary materials.
It is known that there are some spurious connections in the connectivity matrix that should be taken into consideration in the analysis. Also, applying arbitrary thresholding of the FC or SC metrics before calculating the topology properties may influence the result (van Wijk, Stam, & Daffertshofer, 2010). Therefore, we applied a range of thresholds to study the topology properties under different network sparsity. The graph theoretical measurements were calculated over 10 different network sparsity values ranging from 0.05 to 0.5, where the network sparsity is defined as the ratio of the number of edges divided by the maximum possible number of edges in a network. Undirected weighted matrices were used in these calculations.
Four common topological measurements were chosen in this study for analysis of network properties, including the clustering coefficient, shortest path length, global efficiency, and smallworldness.
Detailed information is included in supplementary materials. We evaluated the association between the nodal graph theoretical measurements and the recurrent connection strength w, and the subcortical inputs I, using Pearson's correlation. composite Z-score was the primary cognitive outcome and was created from the linear combination of the Z-scores of the seven cognitive domains measured (executive function, speed of information processing, attention and working memory, learning, memory, verbal fluency, and motor). HAND diagnoses were determined for each participant according to the Frascati criteria (Antinori et al., 2007). The neuropsychological composite Z-score and the seven cognitive domains Z-scores were compared between HIV + BSL, HIV + 12 wk, and HC, respectively. We also investigated the relationship between the graph theoretical measurements with cognitive performance scores.

| Statistical analysis
Comparisons of continuous variables between two independent groups were conducted by two-group Welch's unequal variances t test. Fisher's exact test was used to test any proportional differences in race, gender, and education level between the HIV+ and control groups in Table 1. Pearson's correlation test was used to analyze the association between two continuous variables. A p-value p < .05 was considered statistically significant for a single hypothesis testing problem. For inferential problems that involved multiple hypotheses, Benjamini-Hochberg multiple testing procedure was used to control the false discovery rate (FDR) at <.05 level. The statistical analysis in this study was performed using MATLAB 2017b (https://www. mathworks.com/products/matlab.html) and R version 3.6.1 (https:// www.r-project.org/).

| Demographics
Then, 42 HIV+ subjects were age matched with 46 HC. The HIV+ subjects were scanned before starting the cART treatment, and scanned, on average, after 12-week of the initiation of cART treatment. At baseline, 21 HIV+ individuals had normal cognitive performance, 20 had asymptomatic neurocognitive impairment, and one had mild neurocognitive disorder. The overall cognitive performance, based on the summary Z-score of all cognitive tests, was significantly higher in the HC compared to HIV infected individuals. The mean CD4 count and HIV RNA plasma levels at baseline were 515.8 ± 42.3 cells/mm 3 , and 4.254 ± 0.164 log 10 copies/ml, respectively. After 12 weeks, mean CD4 cell count and HIV RNA levels were 566.4 ± 44.5 cells/ mm 3 and 2.675 ± 2.641 log 10 copies/ml, respectively.

| FC simulation using rMFM
We performed rMFM simulations with 25 different random initializations for each cohort, respectively, 75 simulations in total. Supplemental Figure S1a shows the increase of similarity within 500 iterations for model optimization. The similarities within each cohort are consistent, and all above 0.55, indicating that the rMFM whole-brain dynamic model yields good simulation of FC from SC in the training dataset of each cohort. The similarity between simulated FC and empirical FC shows high consistency across different random initializations.
We found that using different initialization parameters had little effect on the final similarity results (see Supplemental Figure S1b).
The variance of the maximum similarity Z-score for 25 random initializations for HC, HIV + BSL, and HIV + 12 wk are 5.438 × 10 −5 , 8.886 × 10 −5 , and 4.210 × 10 −5 , respectively. Different initialization parameters also had little effect on the final rMFM model parameters (see Supplemental Figure S2), including the recurrent connection w, excitatory subcortical input I, global scaling of the SC G, and neuronal noise σ. Therefore, to minimize this effect, for each cohort, the optimum rMFM model was adopted from the averaged values across the five runs which had the highest similarities.

| Model validation
We validated the rMFM model parameters for each cohort using the test dataset. Then we calculated the simulated neuronal activity for each brain region by feeding the empirical SC from the test dataset into the fitted rMFM model for each cohort. We calculated the simulated BOLD signal for each region by feeding the simulated neuronal activity to the Balloon-Windkessel hemodynamic model. The simulated FC was calculated by Pearson's correlation. We run 1,000 simulations with different random initializations for each cohort using the test dataset. The mean simulated FC matrices for each cohort were calculated as the average across the 1,000 simulations (see Figure 2b

| rMFM model parameters comparison between HIV and HC
The microscale brain dynamic properties of the three groups were investigated by comparing the rMFM model parameters, w and I. Here, we compared these two regional values between cohorts to identify which brain regions changed significantly with a treatment intervention.
The recurrent connection strength w is shown in Figure 3 Figure 4g. We also observed that 9 out of 10 brain regions showed transitions toward HC after 12 weeks of treatment, for example in left and right frontal pole and right precuneus (see Figure 4h). We also observed that the regions which differed in subcortical inputs I were asymmetric (see Figure 4g), with more representation in the right posterior regions. These results are also replicated when using Destrieux atlas (see Supplemental Figure S11). were reproduced using the weighted FC matrix (Supplemental Figure S12). We also investigated the network topology on the empirical SC but did not find any statistically significant difference between groups.

|
Next, we investigated the association between local topological measurements and rMFM model parameters. The nodal clustering coefficient and local efficiency results were correlated with recurrent connection strengths w and subcortical inputs I, as shown in Figure 6.
We found significant correlations between recurrent connection strengths w, nodal clustering coefficient (r = .248, p = .041), and local efficiency (r = .253, p = .037) only in the HC group. These relationships were confirmed using Destrieux atlas (Supplemental Figure S13). 3.6 | Neuropsychological test score comparison, association with FC graph theoretical measurements As shown in Figure 7, HIV + BSL subjects had lower overall composite Z-score on the neurophysiologic test battery, and lower motor function score when compared to HC (uncorrected p < .05  Figure S6).
We also investigated the relationship between Z-scores for different cognitive domains and total Z-score and global graph theoretical measurements. Figure 7c shows the significant linear relationship between neuropsychological test scores and graph theoretical measurements. The correlation coefficients are range from 0.19 to 0.25 with uncorrected p < .05.

| DISCUSSION
In this study, we used rMFM whole-brain dynamic modeling to investigate the local and global brain dynamic changes associated with HIV infection. Our results indicate that HIV infection disrupts microscale brain dynamics in several cortical regions, and HIV antiretroviral treatment improves these brain dynamics.

| rMFM revealed disruption of local dynamics changes due to HIV infection
Among the areas found to have changes in both recurrent connection strength w and subcortical inputs I were the frontal poles, left lingual gyrus, and right lateral occipital areas. The frontal pole regions are part of the prefrontal cortex thus involved in working memory and multitasking (Gilbert et al., 2006). Disruption in these networks has significant behavioral consequences as shown by studies in HIV-infected subjects with active cocaine use disorder (Meade, Lowen, MacLean, Key, & Lukas, 2011). Our study is also consistent with findings that HIV-infection disrupts the frontostriatal circuitry (Heaton et al., 1995).
The lingual gyrus (part of the visual processing system) has also been previously reported to be affected in perinatally HIV-infected youths shown to have reduced regional volume and FC in the HIV-infected population (Kallianpur et al., 2012;McIntosh et al., 2017;Samboju et al., 2018). The balanced integration of excitatory and inhibitory synaptic currents in the local cortex allows the cortical network to operate at high efficiency of information transmission at a low energy cost (Yu, Shen, Wang, & Yu, 2018). Although the rMFM model did not disentangle the contribution of excitatory and inhibitory effect, using rMFM method made it possible to compare local recurrent connection changes between groups.
Among the regions that showed significant differences only in We noticed the maximum similarity Z-scores in HC group were consistently lower than the HIV + BSL and HIV + 12 wk groups, when using training dataset to optimize the rMFM model parameters.
This difference indicates the rMFM model fitting is in favor of HIVinfected group than HC group. We have tried to quantitatively compare the metastability of the three groups using the SD of the Kuramoto order parameter (Acebrón, Bonilla, Vicente, Ritort, & Spigler, 2005;Kuramoto, 1984;Strogatz, 2000). The metastability is not significantly different between groups. We also found the maximum similarity Z-scores are comparable between groups, which indicates the rMFM models can produce a good and comparable representation of three groups.
In our study, we have shown that after treatment, the microscale local brain dynamics improved in both recurrent connection strength w and subcortical inputs I. Specifically, we found 9 out of 10 regions' rMFM metrics in HIV-infected subjects change toward HC levels ( Figures 3 and 4). A similar approach was reported by Deco et al. (Saenger et al., 2017) in the context of deep brain stimulation in Parkinson's disease, but to the best of our knowledge, our study is the first to use rMFM metrics to assess response in brain networks after cART. and traumatic brain injury (Pandit et al., 2013). In our study, we found that reduced smallworldness in HIV + BSL was associated with decreased neuropsychological composite Z-score, executive function and learning Z-score (Figure 7c).
The global network efficiency was also reduced in HIV + BSL, suggesting that the communication between different brain regions We and others did not find significant topology difference in structural network of HIV infected individuals (Samboju et al., 2018;Zhuang et al., 2017). However, two studies (Baker et al., 2017;Bell et al., 2018) have shown disrupted structural networks in HIVinfected subjects compared to HC in terms of reduced global clustering coefficient, global network efficiency, and connection strength.
These two studies differ from our study in the population enrolled.
Our patients had been recently diagnosed and had a relatively high CD4 count at baseline (mean 515.8 cells/mm 3 ) compared to the other two studies. It is likely that persistence of HIV infection overtime causes some irreversible structural CNS changes (Zhu et al., 2013).
The local topological measurements showed some association with the local dynamic properties, specifically, we found the local topological measurements of FC were also correlated with the recurrent connection strengths w and subcortical inputs I (Figure 6 and Supplemental Figure S13). The positive correlations between recurrent connection strength w and clustering coefficient, and local efficiency, respectively, were only found in HC group, suggesting that HIV-infection altered microscale local dynamics, which was also reflected in the change of network topology.
We found the association between microscale local dynamics and empirical network topology using both the Desikan atlas and Destrieux atlas. The HC group consistently shows higher associations than HIV + BSL and HIV + 12 wk (Figure 6a,b, and Supplemental Figure S13a,b). In HIV + 12 wk, a significant correlation was only present in Destrieux atlas. This may be due to the inclusion of more cortical parcels in the Destrieux atlas. The results suggest that the cART treatment tends to improve microscale brain dynamics and the efficiency of the brain information transmission.

| Limitations
There are some limitations of our approach due in part to the available dataset. First, the rMFM modeling requires reliable and accurate SC reconstruction. Our diffusion data have only one shell at b-value = 1,000 s/mm 2 . This single shell scheme may introduce errors for tractography due to fiber-crossing issue (Sotiropoulos et al., 2013). In our study, we have adopted a reliable tractography and SC reconstruction pipeline (Zhang et al., 2018) to mitigate to some extent.
Second, we could only derive a population-based whole-brain dynamic model using averaged FC and SC. The reliability of empirical FC is crucial to the rMFM simulation. A recent study from Patricio group (Donnelly-Kehoe et al., 2019) suggests that in order to derive reliable subject-specific brain dynamics, the total length of rsfMRI acquisition should be 20 min with TR = 2 s. They used both static FC matrices (Pearson correlation of BOLD signal) and dynamic-FC (FCD) metrics to evaluate the simulation results, while in our analysis we only used static FC matrices. Though the FCD estimation requires longer rsfMRI scanning time compared to static FC, our 5 min rsfMRI scan may still limit the static FC reliability to generate subject-specific FC. In this regard, we attempted to derive subject-specific rMFM models using individual SC and FC maps; however, the variance of the similarities was much greater than that of the population-based model. We thought this large intersubject variation for rMFM came from both individual differences and SC/FC noise. It is difficult to disentangle these two effects in our current data. However, we have shown that using group-averaged population-based rMFM modeling is much more stable across different cohorts, and we validated these findings using a different atlas. Thus, the rMFM model parameters (recurrent connection strength w, and subcortical input strength I) and the simulated FC can only be used in comparing groups but not individual parameters, such as CD4 count or cognitive scores. The HIV-infected subjects in our study were mostly male; however, in the HC group the male and female participants were evenly distributed.
Although in the rMFM modeling sex differences between groups were embedded in the final results, we cannot fully assess the possible effect of sex on w and I. Finally, our sample size is relatively small, constraining the stability of rMFM modeling as we needed to split the cohorts into training and test groups. Unfortunately, enrolling cART naïve patients is quite challenging. Future studies using more advanced diffusion scheme and longer rsfMRI acquisition, should provide more reliable subject-specific rMFM model reconstruction.

| Conclusion
We investigated the effect of HIV infection on the microscale local dynamics derived by whole-brain computational modeling. We have identified several brain regions where recurrent connection strengths and subcortical inputs differ among HC, HIV+ untreated, and HIV+ after treatment. Treatment improved local brain dynamics, and this was also reflected in improved brain network topology and cognitive performance. However, short-term treatment did not fully reverse CNS injury. Whether further diverge on CNS injury and cognitive performance occurs between HIV+ and HIV− over longer periods of time will need to be addressed in future studies. In this regard, whole-brain dynamic modeling is a promising approach for assessing CNS injury progression and response to interventions.

ACKNOWLEDGMENTS
The authors thank Alicia Tyrell for providing data management; and Arun Venkataraman and Kyle Murray for the insightful discussion and comments during preparation of this paper. This work was supported by the National Institutes of Health (grant numbers R01 MH099921, R01 AG054328, R01MH118020, and UL1TR002001).

CONFLICT OF INTERESTS
The authors declare no conflict of interests.

DATA AVAILABILITY STATEMENT
The datasets from the current study are available from the corresponding author on reasonable request. The codes are available at https://github.com/yzhuang4/whole_brain_modeling. The codes for rMFM model used in the paper are adopted from https://github. com/ThomasYeoLab/CBIG/tree/master/stable_projects/fMRI_dynamics/ Wang2018_MFMem.

ETHICS STATEMENT
The study was reviewed and approved by the Institutional Review Board at the University of Rochester Medical Center and all subjects signed a written informed consent prior to undergoing study procedures. All subjects underwent a comprehensive clinical, laboratory (chemistry, hematology, and urine analysis), neurocognitive, and neuroimaging evaluation. HIV-infected individuals were assessed before and 12 weeks after starting cART, while HIV-uninfected controls were assessed only once at baseline. Subjects' demographics are listed in