The evolution of disease or the progress of recovery of a patient is a complex process, which depends on many factors. A quantitative description of this process in real-time by a single, clinically measurable parameter (biomarker) would be helpful for early, informed and targeted treatment. Organ transplantation is an eminent case in which the evolution of the post-operative clinical condition is highly dependent on the individual case. The quality of management and monitoring of patients after kidney transplant often determines the long-term outcome of the graft. Using NMR spectra of blood samples, taken at different time points from just before to a week after surgery, we have shown that a biomarker can be found that quantitatively monitors the evolution of a clinical condition. We demonstrate that this is possible if the dynamics of the process is considered explicitly: the biomarker is defined and determined as an optimal reaction coordinate that provides a quantitatively accurate description of the stochastic recovery dynamics. The method, originally developed for the analysis of protein folding dynamics, is rigorous, robust and general, i. e. , it can be applied in principle to analyze any type of biological dynamics. Such predictive biomarkers will promote improvement of long-term graft survival after renal transplantation, and have potentially unlimited applications as diagnostic tools.
Here we present a general framework that allows us to determine such an optimal coordinate or biomarker from longitudinal cohort studies directly, without constructing the Markov state model. The method was originally developed to describe complex dynamics of protein folding [7]–[9]. Briefly, a putative functional form of the reaction coordinate is assumed, for example, a linear combination of features (here the metabolome NMR spectra) that could describe the process. The approach is invariant to the choice of the functional form and the set of observables, provided they contain all the essential information about the dynamics of the process. The coordinate is optimized (trained) on a sample of trajectories, i. e. , realizations of a complex multidimensional dynamical process. This is achieved by choosing the coordinate (e. g. , the coefficients of the linear combination) such that the cut based free energy profile associated with the coordinate is the highest [7], [10]. Namely, given an ensemble of N trajectories () and a reaction coordinate functional form an ensemble of reaction coordinate trajectories is constructed by projecting the multidimensional trajectories onto the reaction coordinate (). The optimal reaction coordinate is found by optimizing the parameters so that the cut-based free energy profile [7], [8] is maximal. , where the partition function equals half number of transitions (crossings) by the reaction coordinate time-series through point y; here and below we set. The CFEP () unlike the conventional histogram based free energy profile () is invariant to reaction coordinate rescaling, insensitive to statistical noise and capable of detecting sub-diffusion. Together they determine the coordinate dependent diffusion coefficient and thus completely specify diffusive dynamics [9]. One can maximize instead the generalized cut based free energy profile where the partition function takes into account each transition through point y with weight equal to the transition distance; for a Gaussian distribution of steps (i. e. , diffusive dynamics) the two optimality criteria are equivalent [10]. If the reaction coordinate is a weighted sum of basis functions, as used here, the optimal values of the parameters (that maximize can be found analytically [10]). In supervised optimization a coordinate that accurately describes the dynamics of transition between two given end states (e. g. , healthy and disease) is determined. Incidentally, the coordinate is the probability of full recovery, i. e. , of ending up in the “healthy” state rather than the “disease” state starting from a current state. It is known as committor or folding probability in protein folding studies [10], [11]. If the two end states are separated by the highest free energy barrier, the transition between them corresponds to the slowest relaxation mode, and an eigenvector, corresponding to the slowest mode is an optimal reaction coordinate [10], [12]. This coordinate can be determined in an unsupervised way without explicit definition of end states.
We have shown that the dynamics of recovery from kidney transplant can be quantitatively described as diffusion on a free energy profile, which is a function of a measurable biomarker. Such a biomarker can be determined in an (un) supervised way from longitudinal cohort studies (patient trajectories), which is optimal in the sense that it is able to discriminate where each patient is on a free energy profile. In particular, the probability of rapid recovery (primary function) can be used to devise optimal treatment. Such an approach is general and can be useful to develop optimal biomarkers for diseases that develop slowly and in a complicated way depending on many factors, or unknown unknowns, such as aging [20], cancer [21], [22] and psychological disorders [23].
Approval was given by the regional ethics approval committee approval number REC Ref: 07/H1306/129 The partition function of the cut-based free energy profile at point equals half the sum of the distances of those trajectory steps that go through point y [14]. More precisely, where (x) is the Heaviside step function and x (iΔt) is the reaction coordinate time series sampled with time interval Δt. The cut free energy profile is defined as and; here we assume that. The optimal coordinate is defined as the one with the highest cut profiles (lowest partition function). The justification of the optimization criteria can be summarized as follows (for the details see the cited references). It can be shown that minimum of, with constrains and, is attained when the reaction coordinate y equals the coordinate - an optimal coordinate [10]. Correspondingly, a sub-optimal coordinate with a lower value of the cut profile, has the mean square displacement which grows slower then linear with time [14]. The latter is an indication that dynamics is not diffusive and that non-Markovain memory effects are at play. Another manifestation of a sub-optimal coordinate, is that it has lower free energy barriers and thus a faster kinetics. The kinetics along the coordinate with the highest cut profile is the slowest [9], [14]. Optimization of the reaction coordinate can be performed in a supervised or unsupervised manner. In supervised optimization a coordinate that accurately describes the dynamics of transition between two given end states (e. g. , healthy and disease) is determined. Incidentally, the optimal coordinate is the probability of full recovery, i. e. , of ending up in the “healthy” state rather than the “disease” state starting from a current state. The optimization is constrained by fixing the value for the coordinate for the two state and [10]. If the reaction coordinate is a weighted sum of basis functions, boundary conditions are given as and where and index the points which belong to A and B states, respectively. The optimal weights () which give constrained maximum to can be found analytically [10]. Here we specified the constrains in the following way, which resulted in a more flexible coordinate. Instead of assuming that each trajectory ends at either 0 or 1, we assumed that each trajectory is constrained to end with either 0 or 1 (in other words a trajectory reaches an end state on the following day). In this case the optimal weights are found by minimizing which equalswhere are indexes that refer to time frame, trajectory and basis function, respectively; the second term of the functional describes the boundary condition with equal 0 or 1 for trajectories connected to 0 or, respectively. The optimal parameters are found by solving the corresponding system of linear equations. To facilitate the visual comparison with the unsupervised results where changed to, which results in a shift and change of scale of the optimal coordinate. In unsupervised optimization the determined coordinate describes the slowest relaxation mode (the second eigenvector) of the stochastic dynamics [10]. If the two states (A and B) are separated by the highest barrier, so the slowest relaxation rate corresponds to the transition dynamics between the states, the second eigenvector reaction coordinate approximates the folding probability (the probability of full recovery here) reaction coordinate in the transition state region - the most important part for the description of the transition dynamics [10], [12]. The eigenvectors can be found by minimizing under constraint [10]. Due to the constraint, the optimization function simplifies to the auto-correlation function. If reaction coordinate is a weighted sum of basis functions, the optimal weights can be found analytically. They are the solution of the generalized eigenvalue problem [10]. The free energy profile that describes the disease dynamics cannot be determined from the patients trajectory simply by computing the cut based (or histogram) free energy profiles because the trajectories are not at equilibrium. The procedure described in Ref [13] was employed. Briefly, assuming diffusive dynamics, the equilibrium free energy profile can be computed from the steady state (non-equilibrium) probability distribution asUsing, and, one obtainswhere and are the cut profiles that measure flux in positive and negative direction, respectively. Note that the method for determining the optimal reaction coordinate was originally derived for equilibrium dynamics; an extension of the framework to non-equilibrium dynamics has been suggested recently [24]. Here we assume that while non-equilibrium sampling affects populations, its main effect on the optimization procedure is in altering the contribution (weight) of the different regions to the optimization functional and can be neglected. The probability of “full recovery” (the folding probability) was computed from the free energy profile asThe success probability with 95% confidence interval were estimated from the trajectories by “add two successes and two failures” approach [25] as, where and and are the numbers of trajectories visiting bin ended up in left or right half of the profile, respectively [26]. NMR spectra were obtained for the water soluble components [27] of erythrocytes taken from 18 kidney transplant patients (up to 9 time points from pre-op to 7 days after surgery). One-dimensional NMR spectra were measured at 499. 97 MHz on a Varian Unity Inova 500 spectrometer at 20°C, using a standard PRESAT pulse sequence. For all samples a relaxation delay of ca. 9 s (three times the longest T1) was applied between scans to allow the spins to fully relax, with 256 transients collected into 16384 data points and a spectral width of 6000 Hz. An exponential line broadening of 0. 5 Hz was applied to each free induction decay (FID) and zero filling to 32768 points was carried out, followed by Fourier transformation. Phase and baseline corrections were carried out using ACD/Labs 12. 0 (Advanced Chemistry Development Inc. , Toronto, Canada) and chemical shifts were referenced to the lactate doublet at 1. 33 ppm. Independently of the NMR data, patients have been divided in three classes based on a clinical assessment of the patients into primary function (PF), delayed graft function (DGF) and acute rejection (AR) with and without primary function. Primary function was defined as immediate recovery of renal function. Delayed graft function was defined as the need for dialysis in the first week following transplantation. Diagnosis of acute rejection was conducted on the basis of biopsy and histological findings. Dialysis was performed on day 5 to patient 4, day 2 to patient 7, day 4 to patient 9, day 7 to patient 17 and day 6 to patient 18. All biopsies were conducted between 6 and 9 days following transplantation. Nine patients had immediate primary function, five patients had delayed graft function and four patients had acute rejection. All but one transplant were eventually successful: in addition to acute rejection, patient 14 also suffered from renal artery stenosis, and the graft was ultimately removed. The immunosuppressive regime and induction agents were the same across the cohort.
The evolution of disease or the progress of recovery of a patient is usually monitored by collecting physical parameters, which may be simply the body temperature for a common cold or properties of tissue samples for e. g. , cancer. Most often clinical decisions are taken based on the current value or because of a sizable change of a relevant parameter. As more advanced diagnostic tools become available, and huge numbers of parameters can be collected at short, frequent time intervals, two related questions arise. The first is, which of the parameters provides relevant information on the progress of disease or recovery as opposed to noise? Is there more information that can be obtained from the history of the evolution of such parameters? Here we propose a novel approach that leads, for the specific case of recovery from kidney transplant, to a positive answer.
physics
biochemistry
research design
proteins
protein folding
protein structure
longitudinal studies
biology and life sciences
physical sciences
biomarkers
biophysics
research and analysis methods
Optimal Reaction Coordinate as a Biomarker for the Dynamics of Recovery from Kidney Transplant