extfieldsofstudy:
Medicine
Computer Science

text:
Tumor Ensemble-Based Modeling and Visualization of Emergent Angiogenic Heterogeneity in Breast Cancer

There is a critical need for new tools to investigate the spatio-temporal heterogeneity and phenotypic alterations that arise in the tumor microenvironment. However, computational investigations of emergent inter- and intra-tumor angiogenic heterogeneity necessitate 3D microvascular data from ‘whole-tumors’ as well as “ensembles” of tumors. Until recently, technical limitations such as 3D imaging capabilities, computational power and cost precluded the incorporation of whole-tumor microvascular data in computational models. Here, we describe a novel computational approach based on multimodality, 3D whole-tumor imaging data acquired from eight orthotopic breast tumor xenografts (i.e. a tumor ‘ensemble’). We assessed the heterogeneous angiogenic landscape from the microvascular to tumor ensemble scale in terms of vascular morphology, emergent hemodynamics and intravascular oxygenation. We demonstrate how the abnormal organization and hemodynamics of the tumor microvasculature give rise to unique microvascular niches within the tumor and contribute to inter- and intra-tumor heterogeneity. These tumor ensemble-based simulations together with unique data visualization approaches establish the foundation of a novel ‘cancer atlas’ for investigators to develop their own in silico systems biology applications. We expect this hybrid image-based modeling framework to be adaptable for the study of other tissues (e.g. brain, heart) and other vasculature-dependent diseases (e.g. stroke, myocardial infarction).

Clinical genotyping has recently demonstrated that subpopulations of cancer cells with unique genomes can exist in different regions or 'microenvironments' of a patient's tumor and that these subpopulations can evolve over time 1,2 . A major driver of inter-and intra-tumor heterogeneity is the underlying heterogeneity of the microvasculature, i.e. the irregular microvascular network and resultant abnormal hemodynamics that also pose a formidable hurdle for developing and delivering effective cancer therapies 3 . Therefore, there is a critical need for new computational tools that enable investigations of tumor vascular heterogeneity and the 'emergent' phenotypic alterations that result from changes in the tumor microenvironment 3,4 . However, in silico investigations of emergent inter-and intra-tumor angiogenic heterogeneity require that we have access to 'whole-tumor' 3D microvascular data as well as such data from 'ensembles' of tumors 3,5 . Until recently, technical limitations such as the lack of high-resolution 3D imaging techniques capable of whole-tumor coverage, insufficient computational power and cost presented major hurdles to such systems biology approaches. Here, we present a hybrid image-based modeling framework to demonstrate the feasibility of recapitulating inter-and intra-tumor heterogeneity of the breast cancer microenvironment in terms of vascular morphology, tumor hemodynamics and intravascular oxygenation.
The origins of microvascular heterogeneity in different organs have been investigated in numerous studies using experimental and computational approaches [6][7][8][9][10][11] . Improvements in imaging resolution now permit acquisition of the 3D microvascular structure and quantification of the morphology of individual microvessel segments 1 within solid tumors [12][13][14] and other tissues 15,16 . Moreover, combinations of imaging and mathematical modeling now permit detailed vessel-by-vessel simulations of microvascular blood flow and oxygen transport 15,17,18 . These advances have been accompanied by the development of new computational approaches for modeling hemodynamics, the transport of oxygen and other molecules such as growth factors and nanoparticle-mediated drug transport in realistic tumor vasculatures 19 . However, traditional image-based computational frameworks are not well suited to modeling the emergent phenotypic heterogeneity that arises from the whole-tumor microvasculature because they are limited by the sparse spatial coverage of the underlying tumor microvascular network (usually of the order of <1 mm 3 ) 15,17,18 .
To overcome these limitations, we developed an image-based modeling framework and specialized data visualization approaches that allowed us to incorporate 3D, high-resolution (~µm) whole-tumor microvascular images in models of inter-and intra-tumor phenotypic heterogeneity of the breast cancer microenvironment. We leveraged a multimodality imaging pipeline that exploits the strengths of complementary imaging modalities such as in vivo MRI, ex vivo micro-CT (μCT) and magnetic resonance microscopy (MRM) 20 to generate an ensemble of high fidelity, high-resolution (~µm) whole-tumor microvascular networks (~GB size dataset) from eight orthotopic breast tumor xenografts. We developed computationally convenient approaches to model hemodynamics and intravascular oxygenation in these whole-tumor microvascular networks. This image-based computational framework made it possible to recapitulate whole-tumor angiogenic heterogeneity in terms of: (i) morphology of the tumor microvascular network; (ii) tumor hemodynamics; and (iii) intravascular oxygenation. We also present novel approaches for visualizing these parameters in 3D space to elucidate the spatial heterogeneity of the angiogenic phenotype within and between whole-tumor microvascular networks. Our tumor ensemble-based analyses of heterogeneity revealed larger dispersion in tumor hemodynamic indices than indices of microvascular morphology. We also discovered that heterogeneity of the angiogenic phenotype led to the emergence of distinct vascular niches within the tumor that uniquely shape its hemodynamic landscape. Finally, we characterized the emergent phenotype arising from these different tumor microenvironments based on their vascular morphology and hemodynamics for the entire tumor ensemble and discovered that hemodynamic parameters alone were insufficient for classifying these different tumor microenvironments. Nevertheless, one could exploit 'ensemble-based' metrics of tumor heterogeneity to classify these microenvironments and identify potential new biomarkers.

Results
Whole-tumor ensemble derived from 3D multimodality imaging data. Figure 1 illustrates a 3D whole-tumor ensemble assembled from multimodality images of orthotopic breast tumor xenografts. The soft tissue contrast acquired from magnetic resonance microscopy (MRM) imaging (~40 µm resolution) is shown in grey and the 3D microvascular network acquired from micro-CT (μCT) imaging (8 µm resolution) is shown in a red color map scaled by the mean vessel diameter. Each 3D whole-tumor microvascular network (~GB size dataset) comprised of tens of thousands of vessel segments and a few hundred thousand nodes. The whole-tumor samples are displayed according to the post-inoculation age they were imaged at (i.e. Day 21 or Day 35) wherein the tumor volumes ranged from 76-450 mm 3 . These unique 3D visualizations provide direct insight into the 3D heterogeneity of the angioarchitectural landscape within and between whole-tumors.
Morphological heterogeneity of whole-tumor microvascular networks. Figure 2 shows 2D, high-resolution maps of vascular morphological heterogeneity in Day 21 and Day 35 tumors in terms of: distance to nearest vessel (D v ), vascular length density (L v ) and vascular surface area density (S v ), respectively. These quantitative maps illustrate the spatial heterogeneity of tumor vascular morphology and help identify unique regions within the tumor. For example, large avascular areas (i.e. areas exhibiting high D v , and low L v , S v ) were visible in the center of a Day 35 tumor (Fig. 2b,e,h), and areas exhibiting low D v (Fig. 2a,b) and high L v and S v (Fig. 2d predictions from image-based hemodynamic modeling match those from conventional modeling approaches. We implemented our blood transport model for a rat mesentery microvascular network with boundary conditions based on experimental in vivo data as described in Pries et al. 21 . We observed excellent agreement between our simulated data and that from Pries et al. 21 , as indicated by R 2 = 0.99 and 0.93 for blood flow rates and discharge hematocrit, respectively (Supplementary Fig. S1a,b). Next, we compared the results of hemodynamic simulations based on our tumor ensemble data with experimental and simulated data from the extant literature. Supplementary Tables S1 and S2 summarize the mean morphological and functional parameters computed for the tumor ensemble along with the corresponding ranges (i.e. mean ± standard deviation) derived from the literature [22][23][24][25][26][27][28][29][30][31][32][33] . The literature data exhibited large variations in morphological and functional parameters due to the varying sizes and tumor types used in these xenograft studies. Wherever possible, we compared our values against data for similar tumor type and size. As shown in Supplementary Tables S1 and S2, the mean values for our morphological and simulated hemodynamic parameters were in good agreement with the ranges reported in the literature.
www.nature.com/scientificreports www.nature.com/scientificreports/ Visualizing the spatial heterogeneity of whole-tumor hemodynamics and intravascular oxygenation. Figure 3 illustrates the spatial heterogeneity of the blood flow rate, Q (ml/s) and intravascular oxygenation PO 2 (mmHg) computed for each whole-tumor microvascular network. We observed well-perfused and well-oxygenated vessels within the tumor rim of Day 21 samples ( Analysis of blood flow rate and intravascular oxygenation heterogeneity in perfused regions of whole-tumors, i.e. regions not devoid of blood flow (Q ≥ 10 −10 µl/s) demonstrated profound inter-and intra-tumor hemodynamic heterogeneity (Fig. 4a,b). Collectively, these data clearly demonstrate the inter-and intra-tumor heterogeneity of perfusion and oxygenation that arises from the heterogeneous microvascular morphologies of these tumors. emergent whole-tumor hemodynamics were more heterogeneous than their underlying microvascular morphology. Figure 5a,b show the coefficient of variation (CV) for each tumor across different vascular parameters. These plots indicate that the distribution of blood flow rate (Fig. 5b) was more heterogeneous than those for the underlying vascular morphological parameters (Fig. 5a). Moreover, the dispersion in blood flow was larger than the dispersion in intravascular oxygen tension by more than two orders of magnitude. Vascular parameters, e.g. diameter D, length L, length density L v , distance to nearest vessel D v and surface area density S v exhibited CV values ≤ 2.5, whereas functional parameters e.g. blood flow rate (Q), shear stress (τ), and velocity (Vel) exhibited CV values ranging from 2-8, i.e. were more heterogeneous. Moreover, we observed

Intra-tumor heterogeneity results in the establishment of unique microvascular 'niches'. Our
simulations demonstrated that the heterogeneity of the tumor microvasculature and resulting hemodynamics combined to establish unique microvascular niches within the tumor. Figure 6 illustrates the four tumor vessel classes identified based on tumor hemodynamic definitions in the extant literature 31,34,35 . Class 1 vessels (Fig. 6a) were hypoperfused and hypoxic (i.e. blood flow velocity, Vel < 50 µm/s and PO 2 < 10 mmHg); Class 2 vessels (Fig. 6c) were hyperperfused and normoxic (i.e. Vel ≥ 50 µm/s and PO 2 ≥ 10 mmHg); Class 3 vessels (Fig. 6d) were either hypoperfused and normoxic (i.e. low velocity and normal oxygenation) or hyperperfused and hypoxic (i.e. normal velocity but low oxygenation). Collectively, Class 3 vessels lack the full functionality of Class 2 vessels, yet are capable of maintaining transport better than Class 1 vessels. Finally Class 4 vessels (Fig. 6g) served as functional 'shunts' (i.e. velocity at least 2× larger than the median blood flow velocity, vessel diameter at least 2× greater than the median diameter, and vascular length smaller than the median length by a factor of 0.5) within the tumor. Figure 6b illustrates how one can visualize this heterogeneous hemodynamic landscape in Tumor 1 by www.nature.com/scientificreports www.nature.com/scientificreports/ mapping each vessel class onto the underlying microvascular network. Almost 90% of the vascular volume in this tumor was accounted for by Class 1 vessels that were poorly or minimally perfused, hypoxic and also exhibited the smallest median diameter among all vascular classes (Fig. 6e). Functionally active vessels of Class 2 exhibited the greatest variation in hematocrit distribution (Fig. 6f), and were localized to the tumor rim along with Class 3 vessels. Class 4 vessels created short pathways of low resistance that diverted blood flow between vessel Classes 1 and 2, and Classes 2 and 3 ( Fig. 6b,g). ensemble-based multiparametric data classify the emergent angiogenic landscape better than either morphological or hemodynamic parameters alone. Finally, we utilized our rich multiparametric, morphological and hemodynamic data to classify the emergent angiogenic phenotype as belonging to D21 (i.e. Day 21 group) or D35 (i.e. Day 35 group). Figure 7a demonstrates that employing a traditional metric such as tumor volume alone was insufficient for classifying the tumor phenotype when hierarchical clustering was used. Figure 7c illustrates how heterogeneity of the hemodynamic parameters precluded such classification. In contrast, Fig. 7b,d illustrate how either the morphological parameters alone, or the morphological and functional variables taken together could classify the emergent angiogenic phenotype more efficiently. It should be noted that one Day 21 tumor (Fig. 1d) exhibited morphological features that overlapped with that of one Day 35 tumor (Fig. 1c), and was thus classified in the cluster corresponding to the Day 35 tumors.

Discussion
Tumor microenvironmental heterogeneity is emerging as one of the most important determinants of drug resistance and therapeutic response in cancer patients 3 . Therefore, to identify more effective drug regimens, a comprehensive characterization of tumor heterogeneity is being sought in the clinic via the collection of multiple biopsy samples from the same patient before and after therapy 1 . While this is a promising start, the search for a non-invasive approach that can be used safely in a large number of patients has resulted in the development of image-based methods for quantifying tumor heterogeneity 36 . Such 'radiomic' approaches include techniques that mine heterogeneity 'features' from large repositories of clinical imaging data for stratifying patients for   www.nature.com/scientificreports www.nature.com/scientificreports/ therapy 37 . Other image-based approaches to investigate the tumor microenvironment include in silico simulations that incorporate imaging data as inputs to mathematical models 38 . The advantage of such computational approaches is that they permit the simultaneous investigation of multiple microenvironmental variables 38,39 . Since access to patient biopsy samples is often limited, high-fidelity imaging data from preclinical xenograft models can play a crucial role in elucidating tumor heterogeneity, especially when used as inputs to image-based computational models of cancer biology. For example, we recently demonstrated the feasibility of using high-resolution, whole-tumor micro-CT data for modeling the blood flow distribution in a whole-tumor microvascular network 40 . However, assessing emergent inter-and intra-tumor microenvironmental heterogeneity requires that we employ 3D data from multiple whole-tumors and interrogate multiple microenvironmental parameters. To achieve this, we developed an image ensemble-based computational and visualization framework that employed data from a collection of tumor xenografts to model the heterogeneity of microenvironmental variables such as vascular morphology, hemodynamics and intravascular oxygenation.
Our tumor-ensemble approach provides a novel method for mapping and quantifying inter-and intra-tumor heterogeneity in whole-tumor microvascular networks. This approach enabled the quantitative visualization of spatial heterogeneity of key vascular morphological parameters at high spatial resolution. We successfully demonstrated how spatial patterns of intra-tumoral heterogeneity could be exploited to identify differences in the vascular phenotype of tumors. For example, we discovered that advanced Day 35 tumors were characterized by large inter-vessel distances, small vascular lengths and surface area densities. These trends were consistent with our current knowledge about tumor progression and the morphological changes that accompany breast cancer angiogenesis 41 . During the highly dysregulated process of tumor angiogenesis, rampant cell proliferation leads to a disproportionate increase in tumor volume relative to functional vessel sprouting, which in turn results in large avascular areas within the tumor such as those observed in our spatial heterogeneity analysis. Moreover, quantifying the distance to nearest vessel provided information on the length scales associated with extravascular diffusion in these irregular microvascular networks. Our 3D distance maps revealed areas of the tumor where effective delivery of drugs and nutrients could prove challenging due to the large extravascular diffusion distances 26 . www.nature.com/scientificreports www.nature.com/scientificreports/ A significant consequence of our image-based computational framework was that it provided a solution to the long-standing technical hurdle of simulating blood flow and intravascular oxygenation in solid tumors using the whole-tumor 3D microvasculature. The profound inter-and intra-tumor heterogeneity of tumor vascular function we observed from our tumor ensemble simulations was consistent with the heterogeneity in blood flow and oxygenation measured in multiple solid tumor samples by other groups 29,31,32,42 . The distributions of blood flow rates and intravascular oxygenation were positively skewed 42,43 , indicating that a large sub-population of tumor vessels was inadequately perfused and under-oxygenated. Perfused vessels within each tumor exhibited high and intermediate oxygen tensions, but the median simulated oxygenation matched microelectrode-based measurements of intravascular oxygen tension in mammary tumors 32,42 as well as other tumors 29,33 (see Supplementary  Tables S1 and S2). While the majority of the simulated blood flow values were within the ranges reported in the literature, we observed that a small fraction (<0.5% of the total vessels) of perfused vessel segments in each tumor exhibited high blood flow velocities (>10 mm/s). This may have resulted from the heterogeneous reconstructed vascular network of tumors or small anomalies in the 3D reconstruction of the microvascular topology from our μCT imaging data. Mapping intravascular oxygenation heterogeneity in whole-tumor networks enabled the visualization of the spatial distribution of perfusion-limited hypoxic vessels, which in combination with diffusion-limited hypoxic regions are known to be major drivers of breast tumor progression and metastasis 44,45 . . Cophenetic correlation coefficients for these clusters were 0.88 0.82, 0.85 and 0.84, respectively. The misclassified D35 sample is Tumor 4 (Fig. 1d) and the misclassified D21 sample is Tumor 3 (Fig. 1c). It should be noted that one Day 35 tumor (Fig. 1d) exhibited morphological features that strongly overlapped with that of one Day 21 tumor (Fig. 1c) and was thus repeatedly found in the cluster corresponding to the Day 21 tumors. www.nature.com/scientificreports www.nature.com/scientificreports/ Moreover, our indices of intra-tumoral hemodynamic heterogeneity could potentially be used to explore how these variables correlate with the probability of metastatic dissemination as well as with treatment outcome 46,47 . A powerful aspect of this image-based computational approach was that our hemodynamic simulations recapitulated emergent properties of the whole-tumor. For example, we were able to demonstrate how heterogeneity of the tumor microvasculature gave rise to unique microvascular niches that comprised of tumor vessels with distinct morphological and hemodynamic phenotypes. Therefore, one can foresee how image-based simulations of these intra-tumoral habitats could be combined with complementary genomic and metabolomics data to better understand the impact of tumor heterogeneity in breast cancer patients 48,49 .
Assessing inter-and intra-tumor heterogeneity using a tumor-ensemble revealed how vascular morphological heterogeneity resulted a heterogeneous functional phenotype. This analysis allowed us to gain critical insights into the emergent angiogenic phenotype. We illustrated how morphological and hemodynamic variables could collectively classify the emergent angiogenic phenotype better than hemodynamic variables alone. Therefore, it might be worth developing a systematic ensemble-based multiparametric heterogeneity index for tumor phenotyping rather than utilizing traditional metrics such as tumor volume 37 . Realistic simulations of angiogenic heterogeneity using the de facto tumor vasculature also have the potential to aide in the identification and development of non-invasive clinical biomarkers of tumor angiogenesis 50 . Collectively, our findings demonstrate the potential of combining structural and functional imaging biomarkers to help characterize the angiogenic status of a tumor or when quantifying its response to antiangiogenic therapies 51,52 .
Despite the advantages of this high-fidelity image-based computational framework, our approach did have some inherent limitations. For example, although we reconstructed more than 90% of vascular network nodes for all tumors to account for incomplete vascular filling of microvessels, we were unable to assess the effect of the remaining 10% of the nodes. Studies using percolation-based network models have shown that a significant fraction (>40%) of vessels needs to be excluded from the microvascular network before any observable impact on network connectivity or simulated blood flow distribution 53 . Since our aim was to model intravascular oxygenation for an entire tumor ensemble, we used a simplified model of oxygen transport 54 . The computational cost of this approach per simulation session was ~10 days on a local workstation (2.8 GHz, quad-core with 8 GB RAM) running parallelized MATLAB code. This simulation run-time could be dramatically reduced by switching to GPU-accelerated computing, an area we are currently exploring.
Since it is challenging to perform 3D reaction-diffusion-convection simulations for the whole-tumor volume, the current formulation for oxygen transport did not include the 3D extravascular diffusion of oxygen in the tissue. Although intratumoral blood flow can be influenced by other factors such as physical stress and the interaction between tumor cells and stromal components 55 , the current computational algorithm did not incorporate this diversity in cellular composition. We have utilized tumor data from two imaging time points (Day 21 and Day 35), however, our computational framework can also be used to study earlier stage tumors or tumors following therapy, as long as each tumor exhibits a patent microvascular network that that can be perfused with a contrast agent that is visible using 3D imaging approaches (e.g. MRM, μCT or multiphoton microscopy). Finally, the current blood flow simulations were conducted for steady-state conditions and therefore do not recapitulate the temporal aspects of hemodynamic heterogeneity known to exist within tumors 14,56 , which is an ongoing area of future development.
It is worth mentioning that our computational approach was based on tumor microvasculature derived from an orthotopic xenograft model of human breast cancer. Since most tumor xenografts only proliferate in immune compromised animals, their microenvironments lack the immune components exhibited by human tumors. Additionally, tumor xenografts do not recapitulate all the genetic diversity observed in primary human tumors 57 . Due to these limitations, our analysis of tumor vascular heterogeneity may not recapitulate every aspect of phenotypic heterogeneity observed in breast cancer patients. However, that was not the primary goal of the current study and our approach did recapitulate salient features of the breast tumor microenvironment such as hypoxia, physiologically anomalous blood flow and the emergence of distinct microvascular niches 3 . We expect future computational applications to be based on images derived from syngeneic tumor models or patient derived xenografts (PDX) in immunocompetent animals. A potential avenue for computationally modeling cancer heterogeneity in patients would be to employ vascular data directly derived from in vivo breast cancer angiography or resected breast cancer tissue in which the vasculature has been labeled by some means. This approach seems feasible owing to recent advances in optical clearing methods in conjunction with novel optical imaging techniques such as light-sheet microscopy. For example, recently a fixed human brain tissue sample (3 × 3 × 1 cm 3 ) was optically cleared to acquire high-resolution images of the brain microvasculature 58 . We are hopeful that these techniques would soon be available for optical clearing of patient derived tumor tissue.
The image-based platform we developed here would be ideal for preclinical investigations of drug development applications, such as studying vascular phenomena at the whole-tumor scale including vessel adaptation and remodeling 56 in response to antiangiogenic therapies, and other treatments that target the vasculature 55 . Multiparametric datasets generated using this platform could also be useful for future clinical investigations. For example, when high-resolution measurements of tumor vasculature in patients become clinically available using in vivo CT angiography 59 , personalized simulations could predict the tumor's angiogenic status and drug distribution levels. Physician-researcher partnerships 49 could facilitate such efforts, since not all clinical settings may have the state-of-the-art computational resources necessary for simulating tumor hemodynamics. As a potential clinical application of this approach, in silico phenotypic data of metastatic breast cancer patients could be combined with their genetic and metabolic profiles to develop patient-specific data-libraries 49 that could be shared widely and be used for predicting metastatic risk or stratifying patient treatment.
In conclusion, we believe that the development of an integrated imaging and computational platform based on high-fidelity tumor ensemble data will help establish: (i) a freely downloadable, multimodality atlas for cancer systems biology investigators and novel in silico applications; (ii) a computational model of antiangiogenic www.nature.com/scientificreports www.nature.com/scientificreports/ resistance that facilitates discovery of translational biomarkers for identifying patients who would benefit most from antiangiogenic agents; (iii) an in silico platform for testing novel therapeutics that can circumvent antiangiogenic resistance in breast cancer. We expect this image-based computational framework to be adaptable to other vasculature-dependent diseases such as stroke 60 and myocardial infarction 61 .

Methods
Imaging of the 3D tumor ensemble. According to the experimental protocol 20 , 3 × 10 6 triple-negative human breast cancer (MDA-MB-231) cells were inoculated orthotopically into the mammary fat pad of eight female Ncr-nu/nu mice. The animals were perfused with Microfil (FlowTech Inc., Carver, MA) a radioopaque vascular filling agent, on Day 21 or 35 post inoculation. Ex vivo μCT imaging was performed at Numira Biosciences (Salt Lake City, UT) on a high-resolution (8 μm) scanner (μCT40, ScanCo Medical, Brüttisellen, CH). The following parameters were used: 55 kVp, 300 ms exposure time, 2000 views and 5 frames per view. Since tumors were perfusion-fixed prior to perfusing with Microfill and imaging was conducted ex vivo 20 , cellular damage from X-ray exposure was not a concern. Ex vivo MRM images were acquired on a vertical bore Bruker 9.4 T spectrometer using a 10 mm volume RF coil (Bruker BioSpin Corp, Billerica,MA) and a T2*w MGE sequence with the following parameters: TE = 4.9/9.7/14. Model for tumor perfusion. To perform hemodynamic simulations, we reconstructed the topology of whole-tumor microvascular networks from 3D µCT images using a bioinformatics-based algorithm we recently reported 40 . The one-dimensional blood flow model 11,40 encompassed the generalized 1D Poiseuille law, mass conservation at vessel junctions and the nonlinear rheological effects of blood flow (i.e. Fahraeus, Fahraeus-Lindqvist and plasma skimming effects) that are known to be significant in the microcirculation 21,62 . An optimization algorithm, that was based in part on a prior study 63 , monitored the vessel segments that feed and drain blood from the network, and optimized the boundary pressures to ensure mass conservation. Using these optimized pressures as boundary conditions, a system of linear equations was solved to obtain pressure distribution at the vessel junctions 40 . For a microvascular network for which complete boundary data is available (e.g. the rat mesentery network used for model validation in this study), our blood flow model can be run without implementing the optimization algorithm.

Model for Intravascular oxygen transport.
We employed the oxygen model introduced by Goldman and Popel with appropriate boundary conditions to estimate intravascular oxygen distribution within each whole tumor 54 . The current implementation is a highly simplified model and does not consider the three-dimensional extravascular diffusion of oxygen in the tissue because it is very challenging to perform 3D reaction-diffusion-convection simulations for the whole-tumor volume. As an approximation, the tumor volume was divided into multiple tissue domains whose volumes are proportional to corresponding vessel segment lengths, L ij . The total flux from the vessel segment is then proportional to tumor oxygen consumption rate; we utilized the rates determined experimentally for breast cancer xenografts by Vaupel et al. 35 . Therefore, this approach is equivalent to a generalized Krogh Cylinder Model, but without the geometric representation; thus, the tissue oxygen distribution cannot be calculated from Krogh-type approximation.
The governing equations for the individual segments can be described as follows: Here, values of blood flow Q b ij , discharge hematocrit H D ij and tube hematocrit H T ij in the individual vessel segments ij of length L ij (where i and j represent branching or boundary nodes) are estimated from the blood flow model. P ij is the oxygen tension in the blood; C bind and S Hb ij are red blood cells' binding capacity and oxygen saturation respectively; J wall ij is the total oxygen flux across the vessel wall for vessel segment of length L ij ; M c is the oxygen consumption rate for breast cancer xenografts derived from Vaupel et al. 35 , m tissue is the tissue mass; n hill is the Hill coefficient for hemoglobin binding cooperativity (2.7) 64 , P 50 is the constant for 50% hemoglobin oxygen saturation (37 mmHg) 64 and α pl is the oxygen solubility in plasma (2.82 × 10 −5 ml O 2 /ml tissue/mmHg) 64 .
www.nature.com/scientificreports www.nature.com/scientificreports/ Equation 1, with application of appropriate boundary conditions, can be transformed to the following system of equations: where, D int and S int are the sparse matrices containing the transport coefficients for dissolved and bound oxygen derived from Equation 1 for branching nodes; P int is the oxygen tension vector for branching nodes; D bound and S bound are the vectors of intravascular transport coefficients for boundary segments; P bound is the oxygen partial pressure vector for boundary segments. This system is solved for oxygen tension (P int ) at all branching nodes and then S Hb values are updated based on Equation 4. The process is iterated until convergence.
Since the oxygen flux per unit length J L ( / ) wall ij ij is assumed to be constant across each vessel wall segment (Equation 2), this term was incorporated only in the calculations for the vessels that transport significant amounts of oxygen subject to following limitations: discharge hematocrit, H D > 0.01 and vascular segment transit time, vstt < 25 s. In a more stringent approach, the local flux per unit length is determined by matching the intravascular and extravascular transport 64 ; nevertheless, the approximation of a constant flux used here provides useful indices of intravascular oxygenation.
Validation of the blood transport model. To validate the predictions from our image-based computational framework, we employed a widely used and publicly available dataset derived from a 546 segment rat mesentery microvascular network 21 . We implemented our blood transport model for this public dataset as described in Pries et al. 65 . The boundary conditions for these simulations consisted of measured blood flow and discharge hematocrit values at all inlets (n = 31), and measured blood flow values at all outlets (n = 5), except the main venular outlet node that was maintained at constant pressure 21 . In addition, we compared the results of hemodynamic simulations based on our tumor ensemble data with experimental and simulated data from the extant literature 21,[27][28][29]31,33 .

3D visualization of heterogeneous tumor vascular morphology, function and vessel classes.
To quantify heterogeneity of the tumor vascular morphology, we generated 3D high-resolution (8 × 8 × 8 µm 3 voxel) maps of distance to nearest blood vessel from the reconstructed whole-tumor microvasculature images in Amira (Visage Imaging, San Diego, CA, USA) using the Distance Map Module. The value assigned to each voxel was its 3D Euclidian distance to nearest blood vessel. We generated 3D maps of vascular length density and surface area density from the reconstructed skeletonized vasculature using custom Matlab code (MathWorks, Inc.). First, every vascular pixel was assigned an effective vascular length (and surface area) based on its position along the vessel. Next, the resulting data was convolved with a 20 × 20 × 20 kernel to estimate average vascular length density, L v (mm/mm 3 ) or surface area density, S v (mm 2 /mm 3 ) in a tumor sub-volume V ROI = 160 × 160 × 160 µm 3 . To render 3D visualizations of spatial distributions of blood flow and intravascular oxygenation, we used the spatial graph visualization module in Amira ® wherein each vessel was scaled by its mean diameter and color-coded by the functional parameter being displayed.
We set thresholds for blood flow velocity (50 µm/s) and oxygenation (10 mmHg) based on tumor hemodynamics definitions in the extant literature 31,35 to identify the first three tumor vessel classes. While vessels exhibiting velocity smaller than 50 µm/s and oxygenation less than 10 mmHg were grouped into Class 1, Class 3 vessels were those that demonstrated greater velocities and higher oxygenation than these thresholds. Class 2 vessels represented the class of vessels that lacked the full functionality of Class 3 vessels, but were still capable of maintaining transport better than Class 1 vessels. Finally, Class 4 vessels represented functional shunts 34 that exhibited velocity at least 2× greater than the median blood flow velocity, vessel diameter at least 2× greater than the median diameter, and vascular length smaller than the median length by a factor of 0.5. To map the spatial distribution of each class in 3D space, we used the spatial graph visualization module in Amira wherein each vessel was scaled by its mean diameter and color-coded by its class.
Calculation of heterogeneity indices. Substantial variations in tumor vascular morphology and hemodynamics make it challenging to compare the heterogeneity of one vascular parameter with another. Therefore, to quantify inter-and intra-tumor heterogeneity within the tumor ensemble, we estimated the coefficient of variation (CV) (i.e. the ratio of the standard deviation to the mean) for each parameter distribution as a normalized measure of its heterogeneity. To determine the strength of the association between the heterogeneity of different parameters, we computed the Pearson correlation coefficient between every parameter pair. statistical Analyses. Two-tailed Mann-Whitney U test (α = 0.05) was performed using the NCSS statistical software (NCSS, Kaysville UT) to determine whether there was a significant difference between the two tumor groups (i.e. Day 21 and Day 35) in terms of median morphological properties such as distance to nearest vessel, vascular length and diameter, as well as median hemodynamic properties such as vascular segment transit time and oxygenation. To ensure valid clustering results, Cophenetic correlation coefficients were computed using the NCSS statistics package (NCSS, Kaysville UT) and values greater than 0.8 deemed significant for hierarchical clustering.

Data Availability
All data supporting the findings of this study are available within the article and its supplementary information files. The Matlab code for all modules of the image-based computational pipeline presented in this study and associated 3D data are available from the corresponding author upon reasonable request.


s2fieldsofstudy: Biology Medicine

text:
Mesenchymal stem cells, implications for pain therapy

History and biology of mesenchymal stem cell (MSC): MSCs were firstly described around 1970s by Friedenstein and his collaborators as a subpopulation of colony forming, clonogenic, plastic adherent stromal bone marrow derived cells, with the capacity to regenerate other tissues (Murphy et al., 2013). Caplan (1991) proposed that these mesenchymal cells were “stem” because of their ability to differentiate to any lineage of mesodermal cells. The richest sources of MSCs are bone marrow and adipose tissue, but they can also be isolated from synovial fluid, periosteum, fetal organs, placenta, umbilical cord blood and amniotic fluid (Murphy et al., 2013; Sherman et al., 2015).

History and biology of mesenchymal stem cell (MSC):
MSCs were firstly described around 1970s by Friedenstein and his collaborators as a subpopulation of colony forming, clonogenic, plastic adherent stromal bone marrow derived cells, with the capacity to regenerate other tissues (Murphy et al., 2013). Caplan (1991) proposed that these mesenchymal cells were "stem" because of their ability to differentiate to any lineage of mesodermal cells. The richest sources of MSCs are bone marrow and adipose tissue, but they can also be isolated from synovial fluid, periosteum, fetal organs, placenta, umbilical cord blood and amniotic fluid (Murphy et al., 2013;Sherman et al., 2015).

MSCs therapeutic use:
MSCs were initially studied and therapeutically proposed only for regenerative aims, due to their differentiation competence, but new correlated functionalities have been recently gaining ground: cell-delivery system and paracrine activity (Brini et al., 2017). Regarding the first function, MSC homing ability is exploited to modulate the release of growth factors and cytokines at the desired site and, moreover, MSCs can be bioengineered by insertion of nucleic acid sequences (e.g., messenger RNA and noncoding RNA), drugs or oncolytic viruses (e.g., myxoma virus to kill glioblastoma cells), to deliver their content to inflammatory areas or tumours, where they can produce or release therapeutic molecules (Sherman et al., 2015). Paracrine activity is linked to the cell-delivery system: the "paracrine hypothesis" states that MSCs are mostly effective through their secretome, which contains trophic, anti-inflammatory, immunomodulatory and antiapoptotic molecules, such as fibroblast growth factor, vascular endothelial growth factor (VEGF)-A, nerve growth factor, transforming growth factor-β and interleukin-10 Chronic pain: Chronic pain is a major public health problem which has detrimental effects on economic and social issues: in the USA, it hits more than 100 million people and it costs more than $635 million per year. Pain is inter-related to biological, psychological and social issues which, in turn, affect quality of patients' life. It often worsens social interactions, reduces leisure activities and limits family contacts. It is bidirectionally linked to low physical activity and to the absence of sleep, which in turn negatively influences the psychological mood of the patient. In the economic framework, chronic pain produces high prolonged costs for the health care system and an extended use of primary * * https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6676889/ 2/4 resources from it; besides, patients' absenteeism from their job may drop the outcome of the company. Neuropathies generally goes together with chronic pain, which in this case is usually treated by antidepressants, antiepileptics and opioids. Opioids misuse and abuse is causing, especially in the USA, drug tolerance and drug-induced hyperalgesia; while the other two lines of medication are limitedly effective. Therefore, there is an urgent need for new treatments and MSCs appeared to be a valid alternative to the existing ones (Dueñas et al., 2016;Hua et al., 2016).

MSC-based treatment and neuropathic pain:
Currently, rationalization of MSC-based therapy is being guided by the "paracrine hypothesis": several studies, carried out on neuropathic pain induced by sciatic nerve ligation, reported that MSC-derived molecules reduce astrocyte reactivity and modulate microglia phenotype switching to an anti-inflammatory M2 behaviour; other papers on diabetic neuropathy proved that MSC secretome is able to significantly lessen neuroinflammation and to restore inflammatory cytokines and signaling molecules balance (Brini et  . We reported that VEGF-A increased in plasma of oxaliplatin-treated rats and rASCs injections resulted in a strong reduction of its plasmatic levels; we also observed that rASC injections brought down pan-VEGF-A and isoform VEGF b levels, which were altered in plasma and spinal cord of treated animals, respectively. Moreover, the isoform VEGF b itself, when injected peripherally, showed a pro-algesic effect, while antibody-mediated blockage of VEGF-A or VEGF b produced an opposite outcome. Our data confirmed the use of MSCs as a novel mechanism of modulation of surrounding cells activity and proposed VEGF-A signaling as a new target for pain relief (Figure 1). VEGF-A mechanism in chemotherapy-induced neuropathic pain modulation.

Conclusions -pros and cons of MSC therapeutic use:
Cellular therapy by MSCs shows many positive aspects: immunosuppressive action which reduces allogeneic lymphocyte response in case of transplantation; minimal ethical concerns; a wide number of tissues to collect MSCs from; an easy in vitro expansion; the ability to autonomously reach the sites of inflammation, due to the expression of chemotactic receptors, where they stimulate tissue repair, release therapeutic molecules or modulate the release of molecules from surrounding cells. Next to the advantages which these putative therapeutic agents have, there are also drawbacks which should be carefully debated and confounding factors which need to be elucidated. First, technical aspects of the experimental studies on this therapeutic tool: the lack of MSCs tracking techniques in humans and, as a result, the consequent difficulties in understanding the fate of the exogenous stem cells and the consequent toxic potential; the absence of common agreement on the isolation method, as many laboratories have their own home-made protocols; different methods of stem cells injection; the correlation of MSC source to a single type of disease or to any of them; the differentiation stage, e.g., a multipotent cell will live longer than a partly differentiated, which will rapidly undergo senescence. Second, issues referring to the therapeutic application itself: different sources of MSCs may correspond to different capabilities; in vitro expansion may lead to accumulation of genetic mutations which could result in cancerogenesis, even if there is no reported case of transformation occurred in vitro (Sherman et al., 2015). In conclusion, this promising therapeutic approach still needs some work around it, to answer these unsolved issues and to ameliorate its use, to exploit it in the most profitable way.
Additional file: Open peer review report 1.

OPEN PEER REVIEW REPORT 1
Click here to view.


text:
Monitoring of miR-181a-5p and miR-155-5p Plasmatic Expression as Prognostic Biomarkers for Acute and Subclinical Rejection in de novo Adult Liver Transplant Recipients

Background and Aims: News strategies for the accurate assessment of the state of immunosuppression (IS) in liver transplant recipients are needed to prevent rejection and minimize drug-related side effects. miRNAs can potentially be used as diagnostic or prognostic biomarkers in transplant patients. This study evaluated the capacity of a plasmatic miRNA panel (miR-155-5p, miR-122-5p, miR-181a-5p, and miR148-3p) as an early non-invasive prognostic and diagnostic biomarker for T cell-mediated acute rejection (TCMAR) and subclinical rejection (SCR) in adult liver recipients. Methods: A total of 145 liver recipients were included. All patients received a calcineurin inhibitor with or without mycophenolate mofetil and methylprednisolone. Plasmatic miRNA expression was assessed by qPCR before and at different time-points after liver transplantation. Results: Seventeen patients experienced TCMAR, and eight were diagnosed with SCR during the protocol biopsy at the 3rd month post-transplantation. Pre-transplantation, miR-155-5p expression was significantly higher in TCMAR patients and in SCR patients than in non-rejectors, and miR-181a-5p expression was also significantly higher in SCR patients than in non-rejectors. Post-transplantation, before transaminase-level modification, significantly increased miR-181a-5p, miR-155-5p, and miR-122-5p expression was observed in TCMAR and SCR patients. Binary logistic regression analyses showed, post-transplantation, that TCMAR risk was better predicted by individual expression of miR-181a-5p (LOGIT = −6.35 + 3.87*miR-181a-5p), and SCR risk was better predicted by the combination of miR-181a-5p and miR-155-5p expression (LOGIT = −5.18 + 2.27*miR-181a-5p+1.74*miR-155-5p). Conclusions: Pre-transplantation plasmatic miR-155-5p expression may be useful for stratifying low-immunologic-risk patients, and post-transplantation miR-181a-5p and miR-155-5p may be candidates for inclusion in early, non-invasive prognostic biomarker panels to prevent TCMAR or SCR and better identify patient candidates for IS minimization. Large prospective randomized multicenter trials are needed to refine the cut-off values and algorithms and validate the clinical usefulness of these biomarkers.

INTRODUCTION
Despite great advances in immunosuppressive (IS) therapy, T cell-mediated rejection (TCMAR), which is the main form of rejection in liver transplantation (LT), has an incidence estimated ∼21-27% (1,2). Moreover, the diagnosis of TCMAR is based on the findings of liver biopsy, which is an invasive procedure that can result in certain complications. On the other hand, the side effects of IS therapy, such as increased risk of renal failure and neoplasia or the impairment of cardiovascular function, has led to strategies for IS reduction to minimize these deleterious events. Accurate methods for stratifying the net state of immune suppression are lacking. The availability of non-invasive biomarkers of an alloimmune response may provide physicians useful information to better identify patients at high risk of rejection and those who are good candidates for IS minimization post-transplantation and to achieve more personalized IS therapy. Monitoring these biomarkers may decrease the requirement for the use of biopsies (3). Furthermore, subclinical inflammatory lesions in the graft, such as portal inflammation with variable grades of fibrosis as described by Londoño et al. or histological signs of rejection without biochemical abnormalities (subclinical rejection; SCR), are less evaluated entities with uncertain significance and prognosis (4)(5)(6).
Several studies have shown the potential role of miRNA analysis as a non-invasive prognostic and diagnostic biomarker of TCMAR (7). The control of gene expression by miRNAs regulates many cellular functions, including differentiation, proliferation, cellular development, and functional regulation of the immune system, and miRNAs also participate in the cellular mechanisms of inflammatory responses (8,9). There is evidence that miRNAs play an important role in post-transplant clinical events, including rejection, disease recurrence, and the development of tumors.
(AR). Serum levels of these miRNAs are increased in patients with liver injury and positively correlate with aminotransferase levels, but this increase occurs before aminotransferase levels are modified during AR, showing the prognostic capacity of these biomarkers (13). In fact, miR-122-5p is the most abundant liverderived miRNA, constituting 70% of the total miRNA in the liver (14), and thus, a significant increase in this miRNA could be associated with hepatocyte damage, toxicity, or viral infection (15). Schmuck et al. also demonstrated that the levels of miR-122-5p, miR-133a, miR-148a-3p, and miR-194 are significantly higher in the bile of liver recipients who develop AR within the first 6 months after transplantation and during the AR episode (16). Recently, Shaked et al., by using the framework of the Immune Tolerance Network Immunosuppression Withdrawal (ITN030ST) and Clinical Trials in Organ Transplantation (CTOT-03) studies, identified in plasma two miRNAs (miR-483-3p and miR-885-5p) that, when combined together in a signature, could diagnose and predict liver rejection with high accuracy [area under the curve (AUC) = 89.5%; 95% CI = 82-96%; sensitivity = 83.8%; specificity = 87.1%; positive predictive value (PPV) = 72% and negative predictive value (NPV) = 93%] (7).
Local inflammation associated with rejection is tightly regulated by the different T helper (Th) subsets, such as Th1, Th2, and Th17. The differentiation of CD4+T cells into these different subsets is strongly influenced by cytokine signaling and by the expression of subset-specific transcription factors, such as miRNAs. Recently, it has been described that miR-155-5p controls the differentiation of CD4+T cells into Th cells (17,18) and participates in the development of regulatory T cells (19). It has been reported that miR-155-5p can also act as a regulator of IFN-γ production in human T and NK cells (17,20). Previous results from multicenter studies have shown that, in peripheral blood, intracellular T cell (CD4 + and CD8 + T cells) expression of IFN-γ correlates with graft rejection and clinical outcome and could be a useful prognostic and diagnostic biomarker in kidney and liver transplantation (21,22). Several groups have demonstrated that miR-155-5p is overexpressed in renal allografts during acute cellular rejection (10,12). Recently, Blaya et al. (23) showed that miR-155-5p expression is altered in both liver tissue and circulating inflammatory cells (peripheral blood mononuclear cells) during liver injury, thus regulating inflammatory cell recruitment and liver damage. Another miRNA closely related to the function of T lymphocytes is miR-181a. Li et al. demonstrated that T cell receptor (TCR) sensitivity and signaling strength can be modulated at the post-transcriptional level by miR-181a. The recognition of donor antigens by recipient T cells in secondary lymphoid organs initiates the adaptive inflammatory immune response leading to the rejection of allogeneic transplants. Allospecific T cells become activated through the interaction of their TCRs with an intact allogeneic major histocompatibility complex; the modulation of selection also argues that this miRNA might directly impact the mature T cell repertoire, which might further affect the onset and/or progression of the T cell alloresponse. miR-181a represents a novel class of regulatory molecules that can modulate the T cell response. Therefore, changes in its expression can regulate the alloresponse against the implanted graft and, consequently, may play a role in the development of rejection (24,25). No data have been reported regarding the potential utility of evaluating plasmatic miRNA expression for predicting the risk of SCR in LT recipients.
The aim of this study was to evaluate the capacity of a plasmatic miRNA panel (miR-155-5p, miR-122-5p, miR-181a-5p, and miR148-3p) as an early non-invasive prognostic and diagnostic biomarker for TCMAR and (for the first time) for SCR in adult LT recipients.

Study Design and Patients
We conducted a prospective, observational study in a cohort of 145 adult patients undergoing LT in our center (Hospital Clínic Barcelona). Combined liver-kidney recipients were excluded. Only patients followed for at least 3 months after LT with plasma samples collected were considered. The study was approved by the Ethics Committees of the center, and all patients provided their written informed consent.

Patient Follow Up and Immunosuppressant Regimens
All LT recipients were managed by transplant hepatologists according to standardized protocols throughout the follow up. After discharge after transplantation, the patients visited the outpatient clinic monthly for the first 3 months and every 2 months thereafter during the first year. Clinical, demographical and laboratory data of the patients were collected. IS regimens were defined by the pre-transplant liver status according to the Child-Pugh classification. Child-Pugh A patients were given a double IS therapy consisting of tacrolimus (TAC) or cyclosporine (CsA) with target trough levels of 8-10 or 150-300 ng/ml, respectively, plus a tapering dose of corticosteroids to be withdrawn at 6 months after LT. Child-Pugh B and C patients, those transplanted due to acute liver failure and retransplanted patients who received induction therapy with a single dose of basiliximab immediately after LT, delayed the start (5th day) of TAC or CsA with target trough levels of 5-8 or 75-150 ng/ml, respectively, mycophenolate mofetil (MMF) 2,000 mg daily, and a tapering dose of corticosteroids to be withdrawn at 6 months after LT. A low dose of corticosteroids (prednisone 2.5-5 mg) was maintained long-term in those patients transplanted due to an autoimmune disease.

Liver Biopsies and Rejection Episodes
Liver biopsies were performed if clinically indicated and if the following criteria were met: aspartate aminotransferase, alanine aminotransferase, alkaline phosphatase or direct bilirubin serum levels higher than 2-fold the upper limit of normal and absence of pathological findings in the abdominal ultrasound examination that could explain these biochemical abnormalities. All rejection episodes were biopsy proven, and severity was defined by the Rejection Activity Index (RAI) Score (26,27). Moderate (RAI 5-6) and severe (RAI > 6) episodes were treated with steroid boluses (500 mg of methylprednisolone for three consecutive days), while mild episodes (RAI 3-4) were treated by increasing the dose and levels of current immunosuppressants. In addition, a per-protocol biopsy was performed at month 3 after LT in all patients in whom a clinically indicated biopsy had not been performed, unless there was a clinical contraindication (mainly, biliary complications). All these patients submitted to a per-protocol biopsy and had normal liver function tests and abdominal ultrasound examination. If histological signs of TCMAR were found, patients were diagnosed with subclinical rejection (SCR). No additional treatment was given for this event, but IS doses were maintained without minimization for at least 4 weeks after biopsy.

Pharmacokinetic Monitoring
Therapeutic drug monitoring of CNi was performed in the Laboratory of Pharmacology. For each CNi, the trough concentration at the 1st week, on the 15th day, and at the 1st, 2nd, 3rd, 6th, 9th, and 12th month post-transplantation were analyzed. Whole blood TAC concentrations were determined by Tacrolimus-CMIA-Architect from Abbot (Wiesbaden, Germany), and whole blood CsA concentrations were measured by Cyclosporin-ACMIA-Dimension from Siemens (Siemens Healthcare Diagnostic, Deerfield, IL, USA) following the manufacturer's instructions. Fresh samples, without having been previously frozen, were analyzed daily. LGC Standard Proficiency Testing was ensured by the participation of our laboratory in the United Kingdom External Analytical Quality Assessment Service.

Plasma miRNA Analysis
At the same time as the clinical visits and pharmacokinetic profiles, plasma miR-155-5p, miR-122-5p, miR-181a-5p, and miR-148-3p expression was assessed by quantitative real-time PCR (qPCR). We selected these miRNAs based on the results of previous studies that showed them as promising biomarkers in the diagnosis and prognosis of rejection and for their involvement in the response mechanisms that the immune system performs against the graft. A total of 4,785 determinations were performed: 1,305 determinations for each miR evaluated (145 patients, 9 visits, and 3 miRs), except for miR-148-3p, which had 870 determinations, because the follow up was performed until 3 months post-transplantation (145 patients, 6 visits, 1 miR). Blood samples (3 ml) were collected into EDTA-K3 tubes at the pre-transplantation visit and at each post-transplantation visit according to the study design. Blood samples were obtained prior to the immunosuppressant administration (pre-dose); at those points concurrent with rejection episodes, the samples were collected before any treatment change was made. After centrifugation (within 2 h) at 3,000 rpm for 10 min, plasma was collected and stored in RNase-free tubes at −70 • C for batched analysis.
Total RNA was purified from 200 µL of patient plasma according to the manufacturer's instructions [miRCURY TM RNA Isolation Kits -Biofluids ref. #300112 from EXIQON (Denmark)]. Briefly, plasma components were lysed with the provided Lysis Solution, and proteins were precipitated with the provided Protein Precipitation Solution. Isopropanol was added to the collected supernatant, and the solution was loaded into the column. The solution was washed with Wash Solutions 1 and 2, the RNA was eluted with RNase-free water, and the concentration and quantity of the total RNA were measured at 260 nm and 280 nm (A260/A280) using a NanoDrop device (NanoDrop Technologies). Hemolysis contamination was tested using a spectrophotometer and measuring the oxyhemoglobin absorbance at 414 nm. An OD scan was performed from ∼200 to 700 nm, and distinguishing absorbance peaks at 414 nm were used to disqualify hemolysis samples. . To monitor the cDNA synthesis reaction, the synthetic spike-in UniSP6 was used for signs of inhibition (prior to the reverse transcription reaction, we added 1 µl of synthetic spikein (108 copies/µl) per 20 ng sample RNA). PCRs were performed using a Light Cycler 480 instrument. The amplification profile was denatured at 95 • C for 10 min followed by 45 amplification cycles of 95 • C for 10 s and 60 • C for 1 min. At the end of the PCR cycles, melting curve analyses were performed. Negative controls with 1 µg of MS2 carrier RNA as a mock template from the reverse transcription reaction were produced and profiled similarly to the samples. The amplification curves were analyzed using Roche LC Software for determination of Cq by the second derivative method. Average Cq values were normalized to the stably expressed reference miR-103a-3p and miR-191-5p, following the manufacturer's instructions. First, the Cq values for all samples were determined; the average Cq of miR-191-5p + miR-103-3p was calculated, and the Cq was calculated as the difference in Cq values between the miRNA target and the reference control (Cq average of miR-191-5p+miR-103-5p). Relative expression levels of target miRNAs were then evaluated within a sample according to the formula 2∧(-Cq), where high values corresponded to higher expression.

Statistical Analyses
Demographic data and results of the molecular analyses were collected in a unified database. Samples were adjusted to a nonparametric distribution. Statistical differences between groups were assessed with the Mann-Whitney test and Kruskal-Wallis test, and correlations between miRNA expression and clinical events were assessed with Spearman's rho test. We used a receiver operating characteristic (ROC) curve analysis to define the optimal cut-offs for differentiating between patient groups with and without TCMAR or SCR. Optimal biomarker cutoff points to discriminate between patients with and without TCMAR were based on ROC curves and calculated with the best Youden index (28) (sensitivity + specificity−1). For the analysis of SCR, we treated the data of this group as an independent group during the evaluation period. Discriminatory capacity was defined by the AUC (0.7-0.8, acceptable; 0.8-0.9, excellent; >0.9, outstanding), with its 95% confidence interval (CI). All analyses were performed using SPSS 23.0 software (SPSS, Inc., Chicago, IL, USA). All data are presented as the median ± standard deviation (SD). A value of P ≤ 0.05 was considered statistically significant. To better evaluate not only the diagnostic capacity of the biomarkers evaluated in this study but also their prognostic utility, we included in the TCMAR box-plot data from patients who exhibited rejection at this time plus the pre-TCMAR data of the patients who had not yet exhibited rejection at this time but who would do so in a later profile. We did not consider data from rejector patients once the TCMAR episode was resolved in the TCMAR box-plot graph. A binary logistic regression model (26) was performed using NONMEM software [version 7.4.1; Icon development Solutions, Ellicott City, MD, USA; (29)] with the Laplacian first order conditional estimation method. TCMAR and SCR occurrence were evaluated as binary data and used as response variables (RVs), with 0 indicating no event, and 1 indicating occurrence of the event. As explanatory variables, miR-155-5p and miR-181a-5p plasmatic expression were used. The probability of the observed score was linked to explanatory variables through the logit transformation to ensure that the estimated probability fell between 0 and 1. Graphic evaluation of the output was performed with R software (30). As a model evaluation, a visual predictive check (vpc) after 1,000 simulations using vpc R package (31) and a bootstrap analysis after 1,000 resamplings using Perl Speaks NONMEM (PSN) were performed (32,33).

Study Patients
From September 2014 to July 2018, 178 patients were included. Twelve patients remained on the LT waiting list at the end of the inclusion period, 6 died before undergoing LT, and 15 patients did not met the minimum follow-up period for several reasons: 5 died before month 3; 1 transplant could not be performed because of a technical impossibility found during surgery; 4 patients withdrew consent; and 5 had no complications but a shorter than 3-month follow up at the time of analysis. The final study cohort consisted of 145 individuals. The main characteristics are shown in Table 1. Most patients were males (72.4%), with a mean age of 56.5 years. The main etiologies of primary liver disease were HCV and alcohol, and hepatocellular carcinoma was the indication for LT in 47.6% of patients. The majority of donors were donors after brain death, with a median

Rejection Episodes
Eighteen episodes of TCMAR were diagnosed in 17 patients (one patient had two events). Regarding the severity, 9 were moderate, and 9 were mild. With respect to the timing, three episodes of TCMAR were diagnosed during the first week posttransplantation, six during the second week, four at the end of month 1, one at month 3, two at month 6, and two at month 12. All of these episodes were recovered with therapy, and no graft was lost due to rejection. Eight patients were diagnosed with SCR at the time of the per-protocol biopsy. Only one of these eight patients developed clinical TCMAR, after the follow-up period expired. As reflected in Table 1, no significant differences were found in the main characteristics of the cohort when comparing those patients who did not develop rejection, patients who presented TCMAR and those with SCR.

Pharmacokinetics
Trough concentrations for TAC and CsA, and MPA and prednisone doses are summarized in Table 2. No statistically significant differences between rejectors (TCMAR or patients with SCR) and non-rejectors were observed in trough concentrations for TAC or CsA neither in MPA or prednisone doses during the period evaluated.

miRNA Expression and TCMAR
Pre-transplantation, no significant differences were observed in the plasmatic expression of miR-122-5p, miR-181a-5p, and miR-148-3p between patients who suffered a TCMAR rejection episode after transplantation from those patients who were free of rejection (Figures 1B-D). Only the pre-transplantation differences between both groups achieved significance for the expression of miR-155-5p (P < 0.001) ( Figure 1A).
Post-transplantation, a significant increase in the expression of miR-155-5p, miR-122-5p, and miR-181a-5p was observed  in rejector patients compared with that in non-rejector patients (Figures 1A-C); however, no significant differences were observed in the plasmatic expression of miR-148-3p ( Figure 1D). Cut-off values for prognostic TCMAR were determined based on the AUC of the ROC curve analysis for each miRNA, with significantly higher levels in patients who showed rejection (TCMAR). ROC curve analysis showed that pre-transplantation, miR-155-5p had an outstanding capacity to discriminate between rejectors and non-rejectors (AUC = 0.921; 95% CI = 0.861-0.981) (Figure 2A). The optimal cut-off value for prognostic TCMAR based on the AUC of the ROC curve for miR-155-5p was 0.075, with 93% sensitivity, 82% specificity, 56.6% PPV, and 100% NPV (Supplementary Table 1).
The results of an analysis of the individual evolution in the expression of miR-155-5p, miR-122-5p, and miR-181a-5p in each rejector patient prior to, during and after the TCMAR episodes showed that the expression of each miRNA progressively increased preceding the TCMAR episode and reached maximum levels at the time of the episode. Once the immunosuppressive treatment was modified to resolve the TCMAR episode, the miRNA levels decreased (Figure 3).
Comparing the miRNA expression evolution with the evolution of serum aminotransferase levels, and taking into account that aminotransferase levels were not normalized until 1-2 weeks post-surgery, in patients who exhibited rejection after 15 days of transplantation (8 of 17 patients), the increase in the plasmatic miR-155-5p, miR-122-5p, and miR-181a-5p expression occurred before aminotransferase levels were modified during TCMAR. Specifically, in patients who showed rejection at 3, 6, or 12 months post-transplantation, these miRNAs levels were above the cut-off value for the risk of rejection in the first weeks post-transplantation, and their increase was earlier and more rapid than the serum levels of aspartate aminotransferase (AST), alanine aminotransferase (ALT), and gamma-glutamyl-transpeptidase (GGT), which usually took place 2-3 weeks before rejection occurred. Furthermore, no significant differences were observed in the evolution of the serum aminotransferase levels between non-rejectors, TCMAR and SCR patients (Supplementary Table 2).

miRNA Expression and SCR
Pre-transplantation differences between non-rejectors and patients who were diagnosed with SCR after transplantation achieved significance for the expression of miR-155-5p (P < 0.001) ( Figure 1A) and miRNA-181a-5p (P = 0.012) ( Figure 1C); however, no significant differences were observed in the plasmatic expression of miR-122-5p and miR-148-3p between SCR patients and those free of rejection (Figures 1B,D).
Post-transplantation, a significant increase in the expression of miR-155-5p, miR-122-5p, and miR-181a-5p was observed in patients with SCR compared with that in non-rejector patients (Figures 1A-C); however, no significant differences were observed in the plasmatic expression of miR-148-3p ( Figure 1D).
In patients diagnosed with SCR during the protocol biopsy, the expression of each miRNA (miR-155-5p, miR-122-5p, and miR-181a-5p) also progressively increased, from pre-transplantation to the 3rd month post-transplantation, and it remained above the cut-off value for the risk of TCMAR during all periods of the study (1 year post-transplantation) ( Figure 5). Instead, in all of these patients, after the 1st week post-transplantation, the serum transaminase levels remained within normal range. One of these patients (Patient #5) finally showed rejection at the 13th month post-transplantation, for which we obtained a sample to evaluate miRNA expression ( Figure 5). While their miRNA levels were already elevated after the immediate post-transplant period, their transaminase values did not increase until 3 weeks before the TCMAR occurred.

miRNA Expression in Patients With SCR and TCMAR
Between rejectors (TCMAR) and patients with SCR, differences in the expression of miR-155-5p were significantly increased in TCMAR until the 2nd month post-transplantation, after FIGURE 3 | Evolution of miR-155-5p, miR-122-5p, and miR-181a-5p expression in patients experiencing TCMAR. An example of the evolution of the miR-155-5p, miR-122-5p, and miR-181a-5p expression in two rejector patient prior to, during and after the TCMAR episodes. The arrow indicates the time of the TCMAR episode.
which no significant differences were observed in either group ( Figure 1A). For miR-122-5p expression, except in the 1st week post-transplantation, no differences were observed in either group ( Figure 1B). Additionally, for miR-181a-5p after the 3rd month post-transplantation, SCR patients showed an increase in expression compared with the TCMAR group, and these differences disappeared 1 year post-transplantation (Figure 1C).

Development and Validation of a Prognostic Algorithm for the Risk of TCMAR and SCR Based on miR-155-5p and miR-181a-5p Expression
A binary logistic regression model was performed to assess the prognostic capacity of the risk of TCMAR or SCR events based on the expression of miR-155-5p and miR-181a-5p individually or in combination as explanatory variables.
The results showed that TCMAR risk was better predicted by the individual expression of miR-181a-5p than by miR-155-5p alone or both miRNAs in combination. Concerning SCR risk prediction, the model that considers a combination of both miRNAs showed a high capacity to prognostic this clinical event. Bootstrap analysis, after 1,000 resamplings, and a vpc, after 1,000 simulations, demonstrated the robustness of the TCMAR and SCR models. The estimated parameters, logit function and bootstrap results of both models are shown in Table 3.
at the end of follow up. We analyzed their miRNA profile expression, and the results showed no significant differences with the remaining non-rejector patients (Figure 6). Individuals with HCV recurrence showed plasma levels of miR-155-5p and miR-181a-5p below the cut-off values established for the risk of TCMAR. In the case of miR-122-5p, the results showed a trend to be higher in the HCV + group with respect to non-rejectors, but this difference was not statistically significant ( Figure 6C). However, when we analyzed the miRNA expression individually, in 3 of 8 HCV + patients the results showed miR-122-5p plasma levels above the cut-off value established for the risk of TCMAR during the period with detectable RNA (Figure 6D). Therefore, for this specific miRNA, the presence of an active HCV infection could be a confounding factor for predicting rejection.

DISCUSSION
To our knowledge, this study is the first to monitor the plasma expression of a panel of miRNAs sequentially in time, pretransplantation and over one year post-transplantation (1st week, day 15th and after the 1st, 2nd, 3rd, 6th, 9th, and 12th month) in de novo liver recipients. This design provides a better evaluation of the changes in the expression of the evaluated miRNAs and its prognostic capacity for the risk of TCMAR. Furthermore, this study is the first to evaluate a panel of plasmatic miRNAs as a potential biomarker for the assessment of the risk of SCR in LT.
This study demonstrates that sequential monitoring of plasmatic miRNA expression levels could reveal useful prognostic and diagnostic biomarkers for TCMAR and SCR during the post-transplantation period. Specifically, we found that the plasmatic expression of miR-181a-5p, miR-155-5p, and miR-122-5p were statistically significant higher in TCMAR and in SCR patients compared with non-rejectors and could diagnose and prognostic TCMAR and SCR with high accuracy. However, no differences in either group of patients were observed for miR-148-3p expression.
Although the results showed that pre-transplantation, miR-155-5p expression was significantly higher in rejectors (TCMAR) and in the SCR group than in non-rejectors, and miR-181a-5p expression was higher in SCR patients than in non-rejectors, and the ROC curve analysis showed an excellent AUC value, the discrete PPV values (∼50%) did not allow for accurately discriminating patients at risk of suffering clinical events from false positive cases. However, the method appears capable of accurately identifying those patients without risk of suffering TCMAR and SCR, because its NPV was 100%, and as a consequence, those patients with an expression below the established cut-off could be candidates to receive IS therapy that is more moderate, minimizing the adverse effects associated with treatment.
Post-transplantation, our results showed a significant increase in the expression of miR-181a-5p, miR-155-5p, and miR-122-5p in patients who experienced rejection (TCMAR). The AUC (>0.940), PPV (>70%), and NPV (>95%) values for these miRNAs were excellent, indicating a very good discriminatory ability to identify patients at high risk of developing TCMAR. Our results are in line with previous reports on the potential of miRNAs as diagnostic biomarkers of AR (10,16,23). In the case of miR-155-5p, our group demonstrated that this miRNA could be a useful prognostic biomarker for the risk of TCMAR in both liver and renal transplant recipients (8). In the case of patients with SCR, the results also showed a significant increase in the expression of miR-155-5p and miR-122-5p, but lower than that in TCMAR patients, compared with the expression in patients free of any clinical event. miR-181a-5p expression in the SCR group reached even higher levels than that in the TCMAR group from the third month post-transplantation; nevertheless, these differences disappeared 1 year post-transplantation. Li et al. demonstrated that the TCR sensitivity and signaling strength can be modulated by this miRNA and thus may be critical for regulating the development of effector cell function (25). The AUC (>0.820), PPV (>70%), and NPV (>90%) values for these miRNAs were excellent, indicating a very good discriminatory ability to identify patients at high risk of developing SCR. The fact that in each TCMAR patient, the individual expression of miR-181a-5p, miR-155-5p, and miR-122-5p gradually was upregulated prior to the TCMAR episode (and earlier than serum transaminases levels were modified during the AR) and reached maximum levels at the time of the episode, clearly shows not only a diagnostic capacity of the miRNAs but also their potential utility as non-invasive prognostic biomarkers. Many previous studies have reported a marked elevation of the transaminase levels in the first hours after liver reperfusion, with a peak on the 1st or 2nd day after transplantation and, unless an episode of acute rejection occurs, they are usually normalized within the first 15 days post-transplantation (34). Therefore, the elevation of these biochemical parameters during the first days post-transplantation is not specific for a rejection diagnosis because there may be other complications that also alter them. It is important to note that, in our study population, those patients who showed rejection exhibited an increase in miR-155-5p, miR-122-5p, and miR-181a-5p expression that was earlier and more rapid than the increase in serum transaminase levels from the first weeks post-transplantation. Furthermore, we observed that in patients with SCR, in contrast to the serum transaminase levels, which remained within the normal range in all the periods evaluated, the expression of miR-181a-5p, miR-155-5p, and miR-122-5p increased progressively from pre-transplantation to the 3rd month post-transplantation. However, in this case, the expression levels kept increasing above the cut-off value established in this study for the risk of TCMAR during the entire period of the study (1 year post-transplantation), indicating that an inflammatory process remained without resolving. In fact, one of these patients ended up experiencing rejection 1 month after finishing the study. Several renal transplantation studies have shown that the detection, usually by surveillance biopsies, of SCR has an impact on the long-term allograft outcome given that SCR can lead to chronic tubulointerstitial damage, chronic graft dysfunction, chronic rejection and graft loss (35)(36)(37). The pathways by which persistent inflammation initiates a fibrogenic response that destroys the renal parenchyma are well-known. In contrast, the role of SCR in the long-term allograft outcome in LT is not clear. Recently, Londoño et al. (5) demonstrated the importance of subclinical inflammatory lesions in long-term stable LT recipients. The clinical significance of these lesions and their response to immunosuppressive treatment is currently uncertain. For this reason, in this context, the clinical utility of a protocol liver biopsy remains controversial, given that it is an invasive method that may favor some severe adverse events, and the benefit that is obtained is not clear; consequently, many centers do not include this procedure in adult LT programs. If this clinical entity could be detected by non-invasive biomarkers, such as the expression of certain miRNAs, the performance of this type of biopsy could be avoided, which would allow for the more accurate selection of patients who are truly candidates for minimizing IS treatment without risk of suffering rejection and for providing the correct treatment of patients with SCR to decrease this alloreactivity.
The logistic regression results demonstrated that miR-181a-5p expression showed a greater prognostic power for the risk of TCMAR than the individual expression of miR-155-5p or both miRNAs in combination. On the other hand, to identify patients at risk of suffering SCR, the model that combined both miRNAs (miR-181a-5p and miR-155-5p) was more robust than the model evaluating the expression of each miRNA. It is evident that a single biomarker may not be sufficient to reflect all of the complexities associated with LT, and it is necessary to identify a panel of distinct non-invasive biomarkers that can detect the degree of alloreactivity of individuals and diagnose and prognostic graft dysfunction as well as injury. According to our results for evaluating the risk of TCMAR, the plasmatic expression of miR-181a-5p (better than that of miR-155-5p) could be a good candidate for inclusion in this panel of biomarkers. To identify patients at the risk of suffering SCR, it would be preferable to include both miRNAs in this panel.
It is important to note that this study evaluated this miRNA panel in an adult LT population that was practically free of HCV infection (only 8 patients who presented HCV RNA positivity at the time of transplantation developed HCV recurrence, and all of them were non-rejector patients, with the remaining patients free of this infection). Therefore, to evaluate the prognostic capacity of these miRNAs for the risk of TCMAR and SCR, clinical confounding factors associated with this clinical event can be discarded. However, it is important to clarify that in HCV + patients, miR-122-5p expression showed a tendency to be higher than that in patients free of rejection, and during the period with active infection, some patients showed levels above the cutoff value established for the risk of TCMAR. Therefore, for this specific miRNA, the presence of an active HCV infection could be a confounding factor to prognostic rejection; in this situation, there is a requirement to test the HCV viral load to rule out patients without risk of rejection. Several previous studies have described significantly elevated serum levels of miR-122-5p in HCV + patients (38,39).
Considering the possible implementation of these biomarkers in the clinical setting routine, some aspects should be addressed. Clearly, the following advantages of plasma markers for the diagnosis and prognosis, in this case, of liver TCMAR or SCR are obvious: less invasive; involves minimal previous manipulation of the sample; and less costly compared with the gold standard method (biopsy). However, there is still more to learn and discover about miRNAs themselves as well as their interaction with their target genes, and, on the other hand, some methodological critical points in the detection of miRNAs should be resolved to avoid discrepancies between centers. There is a need for interlaboratory analytical crossvalidation and comprehensive standardization of all analytical process (e.g., isolation, storage, measurement, and quantification protocols, the selection of housekeeping genes to correct the data to normalize the results, etc.). It is necessary to develop more robust methods and to improve their sensitivity, specificity and reproducibility as well as to define the cut-off values for these promising miRNAs in the context of prospective, randomized, multicenter clinical trials that allow transplant patient stratification associated with a specific clinical event, such as TCMAR or SCR.
Because our analysis is only a single-center study, confirmation by other centers is needed. Furthermore, our study has some limitations. The event size (TCMAR and SCR) is small; in fact, this is a limitation common to all studies of this type of biomarker that has been published to date. It is important to note that the patients in our hospital were treated with triple therapy, and furthermore, this study was the first in which that the study population was practically free of HCV infection. Both of these factors surely impact the allograft response and, consequently, the rate of rejection. However, a large number of samples was included in the present study, which makes us confident in our data and strengthens the importance of our findings because our results still reached significance despite the small sample size. Furthermore, we have internally validated the two algorithms obtained after the binary logistic regression by the bootstrap method, and the results demonstrated the robustness of both algorithms. The study was performed in a Caucasian population, and our findings should also be validated in a separate non-Caucasian population; furthermore, the absence of patients with antibody-mediated rejection (ABMR) in our cohort did not allow us to evaluate the prognostic capacity of these miRNAs for this clinical event. Nevertheless, despite these limitations, this prospective observational study suggests that the sequential pre-and post-transplantation monitoring of the plasma expression of miR-181a-5p, miR-155-5p, and miR-122-5p may well provide key prognostic biomarkers for the risk of TCMAR and (for first time evaluated) for the risk of SCR, and the variation in the level of these miRNAs occurs earlier and more rapidly than that of serum transaminase levels. Pre-transplantation, miR-155-5p expression could be a useful tool to stratify patients at low immunologic risk; furthermore, post-transplantation, the expression of miR-181a-5p and miR-155-5p could be candidates for inclusion in an early, non-invasive prognostic biomarker panel to prevent TCMAR or SCR.

AUTHOR CONTRIBUTIONS
OM performed the miRNA analyses, interpreted the data, performed the statistical analysis, and drafted the manuscript. PR and GC selected and followed the patients. LO and MS performed patient blood samples extractions. PF performed the miRNA analyses. VF performed the pharmacokinetic analysis. LQ performed the binary logistic regression and bootstrap analysis. MN selected and followed the patients and revised the manuscript for important intellectual content. MB conceived and designed the study, interpreted the data, drafted and revised the manuscript. All authors approved the final version of the manuscript for submission.