Atsarina Larasati
Anindya†
a,
Maria-Jose
Garcia-Bonete†
b,
Maja
Jensen
a,
Christian V.
Recktenwald
b,
Maria
Bokarewa I.
c and
Gergely
Katona
*a
aDepartment of Chemistry and Molecular Biology, University of Gothenburg, Box 462, 40530 Gothenburg, Sweden. E-mail: gergely.katona@gu.se; Fax: +46 31 7863910; Tel: +46 31 7863959
bDepartment of Medical Biochemistry, University of Gothenburg, Gothenburg, Sweden
cDepartment of Rheumatology and Inflammation Research, The Sahlgrenska Academy at University of Gothenburg, Gothenburg, Sweden
First published on 14th April 2022
MicroScale Thermophoresis (MST) follows the movement of fluorescent-labelled biomolecules with different sizes along a temperature gradient. The presence of a “contrary trend” pattern, that is, the trend of fluorescence change reversing at higher titrant concentrations, is a well-known problem with uncertain cause. Conventionally, binding curves and kinetic parameters are derived from MST datasets using regression analysis on isolated time windows, while the rest of the data are ignored, and the “contrary trend” fluorescent levels are also usually removed as outliers. This biased approach can be avoided with a more continuous analysis of the entire kinetic process. The Bayesian model of MST progress curves allows the inference of parameters and modelling of the whole experiment. The removal of unusual data points is unnecessary once the anomalous kinetic process is identified. This alternative data analysis approach was applied to our MST datasets from survivin–hSgol2 interactions, and the results show that the binding curves remained sigmoid when all data were included. We were also able to infer the value and uncertainty of the dissociation constant (KD) by ascribing the anomalous data points to a new, linear kinetic component. This approach demonstrates good posterior predictions from the MST process in both short and longer experiments as well as the feasibility of KD inference from short experiments.
An MST experiment is commonly performed by mixing a fluorescent-labelled molecule (target) with an unlabelled molecule (ligand) at multiple concentrations (typically from nanomolar to millimolar concentrations). The mixture is then incubated for some time before it is placed in a capillary tube. Fluorescent labelling is necessary for monitoring the amount of target in the laser focus and observing the time-dependent changes in fluorescence intensity. It is possible to use the intrinsic UV fluorescence of aromatic amino acid residues in the target using a specialized instrument, but that excludes ligands that also produce UV fluorescence (for example, when the ligand is also a protein). Different ligand concentrations would give rise to a wide range of fluorescence intensities and mask the fluorescence of the target.
To create a temperature gradient, an infrared laser is directed towards a region in the capillary tube. As the molecules respond to the temperature change and diffuse away from the spot, the change in fluorescence signal in the illuminated region is recorded over time, resulting in fluorescence time traces. These are normalized to the initial fluorescence intensity (Fnormal). For the determination of the dissociation constant (KD), the target concentration is kept constant and low, while the ligand concentration is varied and acts as the titrant. The dissociation constant is considered to be:
(1) |
This results in a mixture of unbound and bound forms of the target. Selected regions of the raw fluorescence time traces are further analysed, and in an ideal case, the average intensities in these regions result in sigmoidal binding curves, with Fnormal as a function of the ligand concentration.1,3
Currently, few tools are available to analyse the time traces and convert them into a binding curve, including MO.Affinity Analysis from NanoTemper Technologies and the open-source software PALMIST developed by Scheuermann et al.3 Both tools divide the distinctly shaped time trace pattern into separate regions: a steady-state period before the temperature gradient is applied, a short period after the temperature gradient is applied (the “temperature-jump”), and a later thermophoresis steady-state period.3 MO.Affinity Analysis calculates the Fnormal of the binding curve by normalizing fluorescence at the later steady-state to that at the initial steady-state before illumination. The tool also enables advanced users to modify the region range in a time trace dataset. Scheuermann et al. found that calculating Fnormal from (Fh − Fc)/Fc, introduced more noise to their data (where Fh is the fluorescence at the later steady-state or “hot” region, and Fc is the fluorescence at the “cold” early time). They used another approach by calculating Fnormal from Fh/Fc to get a better signal-to-noise ratio in PALMIST.3,4Fnormal is an important parameter because it is the primary data used for modelling the binding curves in the traditional data analysis framework. Although the cold relative fluorescence level is supposed to be close to 1, random and systematic errors affect it. It is important to note that time plays a limited role in traditional analysis; only the fluorescence levels before and after a certain time interval are compared.
Scheuermann et al. later improved the confidence interval for KD estimation with error surface projection (ESP) by comparing it with other methods, such as the analysis of the variance–covariance matrix, Monte Carlo simulation, and bootstrapping. ESP thoroughly checks the error surface considering that defects in a specific parameter may be covered up by other parameters. Users may also choose a suitable combination of the time trace regions for fitting the binding curve according to the expected binding mode of the interacting molecules. The “temperature-jump” region, for example, may not be used to detect molecular interactions that do not directly affect the fluorescent dye.3 Tso et al. added two new models to calculate the Fnormal of molecules with a 1:2 binding mode better, which have been implemented in the most recent version of PALMIST. One model involves splitting the binding curve into two and analysing them separately; this requires more sampling points to cover an adequate range of concentrations in each sub curve. Another model assumes a symmetrical bivalent molecule to account for the cooperative binding sites. Both approaches require the user to have prior knowledge of the expected binding mode.5
While improvements to the analytic tools were made, Scheuermann et al. also documented a peculiar “contrary trend” pattern. In MST experiments carried out to describe the DVD-Actin/VCA′* interaction, the Fnormal value increased as more titrant was added to the sample mix. However, higher titrant concentrations showed a reversal in trend, with the Fnormal decreasing instead. This reversal was also observed in a simulated 1:2 binding interaction when the dataset was fitted to a model assuming a symmetrical bivalent binding partner. Tso et al. noted that trend reversal or a second inflection point at high titrant concentrations might not necessarily come from binding.5 In the case of DVD-Actin/VCA′*, the cause of such a reversal in the pattern is unknown. At the moment, users of both PALMIST and MO.Affinity Analysis may opt to remove the data points manually so as not to skew the fitted curve, which happens when the observations clearly do not follow a random normal distribution.3 Besides introducing a bias, this approach also does not explain why such ‘outliers’ exist, to begin with. An alternative approach is using robust distributions in a Bayesian framework to analyse the MST binding curves, as described previously for survivin and human Shugoshin-like protein (hSgol) interactions.6 For several error models, robust Bayesian curve fitting is visually better and needs lesser repetition to reach a similar level of precision as the standard regression approach. The robust curve fitting of survivin–hSgol interactions shows that the model captures the central tendency of the data even without removing the contrary data points.6
Progress curve analysis is mostly associated with recording the progress of the chemical reaction dynamics and solving equations to describe continuous reaction kinetics. It has the advantages of inferring parameters using the whole kinetic curve and requiring fewer repeats to estimate these parameters. The limited usage of progress curve analysis is often attributed to its complex mathematical modelling. Increasing advancements in computational power, however, have enabled the integration of rate equations and optimization of their parameters to fit the experiment progress curves.7,8 The optimization of parameters and estimation of the uncertainty of parameters can be performed in a Bayesian framework, which is uncommonly used for inference in progress curve analysis. The Bayesian inference has recently been adapted to describe complex biophysical systems; it uses a continuously self-refining model as more data are added to improve the accuracy and precision of the model.9 By learning from a wide variety of training datasets, the model progressively makes better estimations on the actual dataset. The continuous learning nature of Bayesian inference makes it highly suitable for pairing with the diverse experimental conditions in progress curve analysis. Choi et al. showed that, after applying Bayesian inference on enzyme reaction dynamics with diverse kinetics, the bias in the progress curves was significantly reduced compared with conventional regression. The conventional model showed a rise in errors in the posterior samples with increased enzyme concentrations. Meanwhile, the posterior samples from the progress curve model were consistently more accurate.10 The ideal process of thermophoresis is not equivalent to the changes in concentration in a chemical reaction for which progress curve analysis is applied previously. The empirical progress of diffusion and other (photochemical) processes can be modelled with independent exponential (and linear) components even during a thermophoresis experiment. Therefore, we keep using the term progress curve analysis in a general sense and not limited only to modelling chemical reaction kinetics.
Here, we show that a certain type of anomalous MST data can be successfully modelled with the addition of a single linear kinetic parameter together with the exponential process that is typically attributed to thermophoresis and binding-dependent fluorescence change. Through progress curve analysis, we have identified a new kinetic component that corresponds to the contrary trend observed previously. By applying Bayesian progress curve analysis to our MST data, we were able to extract the kinetic parameters without having to remove any outliers, as well as determine more precise KD values than those determined by traditional regression analysis. We could also quantify the uncertainty of KD obtained from the hierarchical Bayesian framework. Our results show that the fitting of MST traces from short experiments had a similar quality to those from longer ones and presented reproducible posterior sample predictions. Although the KD posterior probability distribution of the short experiments was understandably broader, the location of the posterior peak was identical to that obtained in longer experiments. In practice, this would allow us to use shorter MST time traces to estimate KD.
Fast decay is usually associated with a temperature effect on the yield of the fluorophore itself, which may or may not be affected by the altered chemical environment upon complex formation. In the shown example (Fig. 1a), the fluorescent labels are chemically attached to survivin, and the peptide could make direct contact with the fluorophore. The slow decay is assumed to represent the thermophoretic movement, which is a diffusive process. It is modified by the size, charge, and shape changes occurring in the protein and/or in the hydration shell upon complex formation.
Fig. 1 Experimental MST data. (a) Primary thermophoresis data from serial dilutions of hSgol2 peptide when incubated for 5 min with fluorescent chemical-labelled survivin. The cold (−3 to −1 s) and hot (27–29 s) regions used to analyse the thermophoresis binding curves are represented by blue and red vertical lines, respectively. (b) Thermophoresis binding of hSgol2 at varied hSgol2 concentrations. The blue dots represent measurements obtained from a dilution series.6 |
There are many potential advantages of the (Bayesian) thermophoretic progress curve analysis. For example, one can estimate the exponential amplitudes from a continuous curve without limiting the analysis to the comparison of arbitrarily selected regions in the curve (Fig. 1). The early time points are more valuable and carry information about both exponential processes. Thus, one may even consider ending the experiment before a steady-state fluorescence level is reached. A trivial advantage is that one can perform the experiment faster, but it may also be practically impossible to reach a steady-state fluorescence level. Besides those caused by photobleaching, there are other yet unexplained drifts (modelled with linear time dependence in the absence of a better kinetic model) that continue to act on longer timescales after the exponential T-jump and even when the thermophoretic processes are practically over. We show evidence for this by modelling the progress curves in Fig. 2a.
In Fig. 2, we have compared the fit of a short experiment curve (9 s; Fig. 2b) with a fully recorded (30 s; Fig. 2a) progress curve and the inferred KD distributions. As expected, the probability distribution of the exponential amplitudes and KD parameters inferred from the 9 s curves were broader compared with the 30 s curves, but for most practical purposes, 9 s is sufficient. The mean of KD posterior distribution determined for the long experiment was 1.8 μM (s.d. 0.18 μM), which is in line with a robust analysis method described previously on a triplicate of observations (1.6 μM).6
In Fig. 3, we have compared the influence of ligand concentration on the magnitude of the linear kinetic constant in a fully recorded experiment. At 2 μM hSgol2 peptide concentration, an anomalously large positive gradient value was seen, which forces the progress curve to keep increasing without reaching a steady-state level. At 30 s, the Fnormal value was higher than the presumably fully unbound survivin when 15 nM hSgol2 peptide concentration was applied.
Disregarding this linear background process described by one parameter per curve, the total exponential amplitudes were perfectly sigmoid, as required by the theory (Fig. 4), while the posterior predictions remained remarkably good for all progress curves in the experiment (Fig. 2a). The nominal value of Fnormal may look very unusual towards the end of the experiment, but the total amplitude of the exponential components, and their rates were in line with the equivalent parameters at adjacent ligand concentrations. The linear components showed a concentration-dependent tendency: the gradient increased above 0 both towards the high and low ends of the ligand concentration range. The observation of a near-linear component is not rare. Scheuermann et al. showed a “contrary trend” when studying the DVD-Actin/VCA′* interaction.3 In their study, Fnormal returned to an intermediate level in a strong near-linear manner at high VCA concentrations. A common practice is to remove these contrary data points at high ligand concentrations from the binding curves. Unfortunately, this procedure does not solve the problem, because the “contrary trend” is already present at other concentrations. We could see this smooth variation directly in our data and, if left untreated, this linear trend influences the obtained KD values. A routine solution is to make the analysis of binding curves robust against seemingly anomalous occurrences.6 Alternatively, by thermophoretic progress curve analysis, one can attempt to isolate the interesting kinetic components of thermophoresis from the kinetic background. The process of learning is not only to recognize but also to be able to focus on the relevant core details. The hierarchical nature of this Bayesian network makes it especially sensitive to exponential processes and sigmoid concentration dependence of ligand binding. Every other process that does not fit this framework is not considered relevant. Progress curve analysis is powerful but requires more computational resources, and the MST progress curves sometimes defy even very general expectations. Fortunately, in the described case, a linear time-dependent kinetic background could be used, which already extends the potential range of the curve.
Fig. 5 shows the comparison of the exponential phase amplitudes inferred from the MST data involving the interaction of chemically labelled survivin and the hSgol2 peptide. Both amplitudes showed a correlation with the ligand concentration, and for the fast phase, the transition appeared at a lower ligand concentration. The magnitude of rate constants did not seem to correlate with the ligand concentration in both the fast and slow exponential phases.
Clearly, there is room for further improvement of the thermophoretic progress curve analysis. We were only partly correct with our assumption that the linear background process is common in the pre- and post-heating phases of the experiment. In most cases, the preheating trace is not necessary for the accurate estimation of the post-heating trend; however, it may be worth linking them together with a common hyper-parameter.
After serial dilution and incubation for 5 min, the experiments were performed using 20% MST power and 40% LED power. The MST traces were recorded using the standard parameters: 5 s MST power off, 30 s MST power on (these two periods were used for the progress curve analysis) and 5 s MST power off.
K D is a global variable with the same a priori expectations as described previously (eqn (2)).6 Likewise, the probability distribution of fluorophore concentration is identical (eqn (3)). U and B represent the pure total amplitudes of the exponential processes for the unbound and fully bound fluorophores, respectively. It is also possible to change the definition of U and B to link them to the amplitude of only one of the exponential components.
(2) |
(3) |
(4) |
(5) |
The a priori assumptions for U and B were identical and described by a flat β distribution. By assuming positive amplitudes, we ensured that fluorescence would not increase beyond the initial level due to exponential (T-jump and thermophoretic) processes (we leave that possibility open through a linear process though). Not fully bound fluorophores at ligand concentrations cn have intermediate exponential amplitudes Atotal,n determined by the law of mass action, and are entirely the deterministic combination of the stochastic components above:
(6) |
The experimental data were modelled as part of a normal distribution, and its scale parameter ε was a single global variable with a lognormal a priori distribution in our model. This choice was motivated by the belief that the errors do not vary from capillary to capillary.
(7) |
Local random variables were linked to each of the n thermophoretic progress curves, and their models consisted of one linear and two exponential components. The linear process was modelled from the beginning of the fluorescent signal recording, and the exponential processes started from IR laser irradiation:
Ln(t) = ν0,n(t + 5 s) | (8) |
E1,n(t) = A1,ne−ν1,nt | (9) |
E2,n(t) = A2,ne−ν2,nt | (10) |
Since Atotal,n = A1,n + A2,n, they can be linked together with a single, curve-associated random ratio parameter (Rn) varying between 0 and 1. Rn can be conveniently modelled with a β distribution, and a slightly asymmetric prior expectation would ensure that exponential processes with larger and smaller amplitudes are grouped together for comparison.
E1,n(t) = RnAtotal,ne−ν1,nt | (11) |
E2,n(t) = (1 − Rn)Atotal,ne−ν2,nt | (12) |
(13) |
Since the progress curves were normalized to 1 at t = −5.0 s, the final fluorescence level approached by the two exponential components in an ideal experiment (i.e., an experiment without the linear component) was In = 1 − Ln(5.0 s) − Atotal,n.
The remaining rate parameters were modelled as random variables with the following a priori parameters for the uniform and lognormal likelihood functions:
(14) |
(15) |
(16) |
Before t = 0 s, the progress curves were modelled as:
(17) |
When t ≥ 0 s, the progress curves were modelled as:
(18) |
Before sampling the posterior parameter space, the model was scaled by the Automatic Differentiation Variational Inference (ADVI)11 algorithm, and subsequently, 7000 samples were collected in four parallel chains using the No-U-Turn Sampling algorithm (NUTS).12 Notably, ADVI was not the default scaling method, and the default jitter+adapt_diag method failed to provide a suitable starting point. The model parameters in the last 5000 samples exhibited steady-state fluctuations around the a posteriori maximum and converged to identical distributions in the four parallel sampling processes (Fig. 7). For the final analysis, the 5000 samples from the four chains were merged to yield 20000 samples.
The computational speed was around 1.9 iterations per s on a Linux workstation (i7-3970X CPU at 3.50 GHz clock frequency), with four parallel chains growing at each iteration. An interactive implementation is available through Google's Colab, on the link: https://colab.research.google.com/github/Katona-lab/MST-analysis/blob/main/MST_progress_curve_analysis.ipynb.
At a typical load in a CPU-only Colab environment, the MCMC sampling performance was approximately 1.6 iterations per s, and it was growing only one chain at a time. The posterior samples were also used for 1000 new predictions. These posterior predictive samples were demonstrated by their median and their highest density interval (95%) at each time point on the progress curve.
Footnote |
† These authors contributed equally to this work. |
This journal is © The Royal Society of Chemistry 2022 |