Bayesian inference of the viscoelastic properties of a Jeffrey’s fluid using optical tweezers

Bayesian inference is a conscientious statistical method which is successfully used in many branches of physics and engineering. Compared to conventional approaches, it makes highly efficient use of information hidden in a measured quantity by predicting the distribution of future data points based on posterior information. Here we apply this method to determine the stress-relaxation time and the solvent and polymer contributions to the frequency dependent viscosity of a viscoelastic Jeffrey’s fluid by the analysis of the measured trajectory of an optically trapped Brownian particle. When comparing the results to those obtained from the auto-correlation function, mean-squared displacement or the power spectrum, we find Bayesian inference to be much more accurate and less affected by systematic errors.


Scientific Reports
| (2021) 11:2023 | https://doi.org/10.1038/s41598-021-81094-x www.nature.com/scientificreports/ of a particle in a viscoelastic fluid exhibits memory effects, i.e., it is a non-Markovian process 23,24 . Nevertheless, Markovianity can be regained at the cost of additional degrees of freedom describing a multivariate Ornstein-Uhlenbeck process 32,33 . For such stationary, Gaussian and Markovian processes, the likelihood function which corresponds to discretely observed sample paths can be exactly calculated in terms of the process parameters. Using the Bayes' theorem, the probable parameter values concealed in given discrete observations can be obtained with the corresponding probabilities 33,34 . Such a Bayesian theory based time domain method is thorough, accurate and avoids limitations such as non-trivial Fourier transformations, data range dependent uncertainty observed in the measurement from the auto-correlation function (ACF), power spectral density (PSD), mean-squared displacement function (MSD) etc.
Here, we demonstrate a Bayesian inference based technique to evaluate the rheological properties of a Jeffrey's fluid from the partial observation of the position fluctuation of an optically trapped Brownian particle. We calculate the likelihood function related to the probe's trajectory and estimate the most probable parameter values and their standard errors by numerically maximizing the likelihood with respect to the parameters. Further, we compare our method with the measurements from the position ACF and show that the estimations obtained from the proposed method are more reliable compared to the regular measurements from the ACF. Note that PSD and MSD are directly linked to the ACF, and hence for comparison, the analysis of only one among these is sufficient. We also show numerically and experimentally that both the methods lead to the same estimations for a given time series provided a least-squares fitting is conducted carefully over a certain data range. Furthermore, we also investigate how the length of the time series and the sampling time step affects the accuracy of the rheological parameters, similar to recent work 35 where authors quantitatively determine the efficacy of their method.
Finally, we compare our results with macroscopic rheological measurements for various viscoelastic fluids and obtain good agreements.

Theory
We describe the one dimensional Brownian motion of a microscopic particle of mass m confined in a harmonic potential of stiffness k in a homogeneous and isotropic viscoelastic medium by using the generalized Langevin equation (GLE) 24,32,34 : where x(t) is the position of the particle at time t, Ŵ(t − t ′ ) is the friction kernel and η(t) is the stochastic noise which represents the thermal agitations of the fluid molecules. Within the Jeffrey's fluid model (generalized Maxwell model with single relaxation time), the memory kernel takes the form 9,24,32,36 and to satisfy the fluctuation-dissipation theorem (FDT) (and hence causality), the noise correlation function η(t)η(t ′ ) = k B TŴ(t − t ′ ) ; k B is the Boltzmann constant and T is the absolute temperature 32 . The first term in the expression of Eq. (2) is due to the purely viscous solvent and the second term corresponds to the polymer network present in the fluid. The noise η(t) can be represented as η(t) = η 0 (t) + η 1 (t) , where η 0 (t) and η 1 (t) are two independent Gaussian random processes with mean zero and correlations η 0 (t)η 0 (t ′ ) = 2k B Tγ 0 δ(t − t ′ ) and , respectively. Note that the dynamics of the particle at time t depends on its history and thus the process is non-Markovian. At time scales, where inertia is negligible (overdamped limit), Eq. (1) can be written as 32,37 In order to model the system as Markovian, we introduce an auxiliary variable we can write where φ 0 and φ 1 are two independent delta-correlated Gaussian random variables with mean zero and standard deviation one. If we consider, (4) becomes a multivariate Ornstein-Uhlenbeck process 33 and the transition probability density of the state Y = x X , i.e., the probability density P 1|1 (Y ′ , t ′ |Y , t) of a transition from Y at time t to Y ′ at time t ′ , follows a Fokker-Planck equation (FP) ∂ t P 1|1 = L P 1|1 , with the Fokker-Planck operator The solution of the FP is a multivariate normal distribution with mean µ(�t) = exp (− �t)Y and covariance �(�t) = σ − [exp (− �t)σ exp (− T �t)] ; t = |t ′ − t| . In the limit of t → ∞ one obtains the stationary probability distribution P 1 (Y ) which clearly has a zero mean and covariance σ = YY T . , σ and D are related by the stationary condition L P 1 (Y ) = 0 , which results into σ + ( σ ) T = DD T33 . A straightforward calculation yields 32 The auto-regressive (AR) representation of the process, thus, can be written as 34 where ǫ n is a zero mean Gaussian random vector with co-variance (�t).

Bayesian inference I
If we consider the observations of the time series with N sample points Y ≡ (Y 0 , Y 1 , Y 2 , ..., Y N ) being independently selected from the stationary distribution P 1 (Y ) , then according to the Bayes' theorem if is a set of unknown parameters, the conditional probability of the parameters given the observations, P 1 ( |Y ) (posterior probability) is proportional to the product of the conditional probability of the observations given the parameters, P 1 (Y | ) (likelihood function) and the prior probability of the parameters P( ) 16,33,38 , i.e., Therefore, for a non-informative prior probability, ln P 1 ( |Y ) ∝ ln P 1 (Y | ) . Straightforwardly using the Gaussian nature of the stationary probability described above 16,33 , The maximum a posteriori probability (MAP) estimation is obtained by optimizing ln P 1 ( |Y ) with respect to . The corresponding standard error is obtained from the square root of the inverse of the Hessian matrix, H = −∇∇ ln P 1 ( |Y ) . Here, if the parameter is σ then the MAP estimation is Clearly, from Eqns. (7) and (11) we conclude with the corresponding standard error k * √ N . However, this method can not provide all rheological parameters involved in the process.

Bayesian inference II
Because of the Markovian property of our process, the conditional probability of the whole sample path given a set of unknown parameters is 16,33,38  www.nature.com/scientificreports/ and similar to Bayesian inference I, the log posterior probability ( ln P( |Y ) ) for non-informative prior information is proportional to the log likelihood function ( ln P(Y | ) ). Nevertheless, the MAP estimation by the optimization of ln P( |Y ) depend both on the particle's position x and the auxiliary variable X 33 . Therefore, we need an appropriate likelihood function which can be calculated only using the partial observation of the process, i.e., the trajectory of the particle x ≡ (x 0 , x 1 , ..., x N ).
Likelihood function. We can use the Kalman filter method 34,39,40 to calculate the likelihood function from the partial observation of the process Y . For the n-th step, one can write where C = c 0 ; c is the calibration factor related to the experimental time-series which can be set to unity if the data is already in desired units. Using the AR structure (Eq. 8) and the partial observation of the process, the Kalman filter algorithm iteratively approaches to the likelihood P(x 1 , ..., x N |x 0 , ) for a certain set of the unknown parameters . The log-likelihood function has the form 34 where x n|n−1 and n|n−1 are the conditional mean and variance of x n given x 1 , x 2 , ..., x n−1 and . Further, as discussed above, the optimization of the log-likelihood function with respect to the parameters (we are interested in the set of parameters = (k, η 0 , η 1 , τ 1 ) ) can yield the MAP estimations for the parameters. The corresponding standard error can be calculated using the Hessian matrix. For convenience, we use the Bayesian I to infer the MAP estimation of the trap stiffness, i.e., k * (Eq. 12) and further optimize the log-likelihood function (Eq. 14) for = (η 0 , η 1 , τ 1 ) . Because of a large amount of literature discussing the Kalman filter 34,39,40 , we describe it briefly in our context.
Kalman filter. The Kalman filter algorithm for a linear system generally takes an observation model which linearly depends on the state of the system. Here, the state variable follows the linear transition equation where F = exp (− �t) and ǫ n ∼ N(0, �(�t)) ; the observable x n depends linearly on the state variable as in Eq. (13). Further, it brilliantly utilizes the Bayes' theorem to recursively approach the likelihood function P(x 1 , ..., x N |x 0 , ) . According to Bayes' theorem, the probability of Y n given the observation x n , the previous state Y n−1 and the parameter set is where P(x n |Y n , Y n−1 , ) (n-th iteration for the likelihood function) is the probability of the observation x n given Y n , Y n−1 and ; and P(Y n |Y n−1 , �) (n-th iteration of the prior) is the transition probability of the state Y n−1 to Y n in the sampling time step t , given the parameter set , which is a normal distribution as discussed in the previous section. Optimization of P(Y n |x n , Y n−1 , ) should infer the expected n-th step of the state given the observation x n and parameters , i.e., and the corresponding covariance, i.e., Straightforwardly, we can update the prior probability distribution function to another normal distribution of mean Ŷ n+1|n = FŶ n|n and covariance ω n+1|n = Fω n|n F T + � . The updated likelihood probability distribution function, i.e., the probability of x n+1 given the previous observations ( x n , x n−1 , x n−2 ,...x 0 ) and , evidently, is a normal distribution of mean x n+1|n = CŶ n+1|n and variance n+1|n = Cω n+1|n C T . Therefore, the implemented Kalman filter algorithm is the following: Therefore, to implement this method one should use the iterations Eqns. (17) to calculate the log-likelihood function Eq. (14) and optimize with respect to the desired set of parameters.

Numerical analysis and comparison
We numerically solve Eq. (5) with input parameters τ in 1 = 1 s, η in 0 = 0.001 Pa s, η in 1 = 0.1 Pa s, k in = 0.1 × 10 −6 N/m and analyze the simulated data to study the accuracy of the method in comparison to the fitting methods. Note that the smallest time scale of the process τ S = 1/ k γ 0 + γ 1 γ 0 τ 1 and the largest is τ 1 . The input parameters are chosen so that they are experimentally feasible and τ S becomes 0.01 s for convenience. In Fig. 1a-d, we show the estimations of the parameters for several sampling time steps (normalized by τ S ) with varying total sampling time ( T s ) normalized by τ in 1 . We estimate the trap stiffness using the Bayesian I method and use that in Bayesian II to evaluate other parameters. Note that for �t/τ S ≤ 1 , the MAP estimations of the parameters reach their true values with ∼ 1% standard deviation when the normalized sampling time T s /τ 1 crosses ∼ 500 . Further increase in T s /τ 1 , indeed, improves the precision of the result. However, if �t/τ S > 1 , i.e, if the sampling time step is lower than the smallest time scale of the process, T s /τ 1 should be greater than ∼ 2000 for reliable estimations.
Auto-correlation function. The position auto-correlation function (ACF) of a Brownian harmonic oscillator in a viscoelastic Jeffrey's fluid, as described in the theory section, is given by 24,32 Ideally, the measured ACF from the position time series converges to the functional form (Eq. 19) in the limit T s → ∞ . However, in practice, we can assume that the calculated ACF for each time lag ( t i ) is Gaussian distributed around the trend f (t i ) (Eq. 19) with variance S, and all the measured ACF data for different time lags are independent of each other. Thus, the corresponding likelihood function is where N(f i ( ), S) represents a normal distribution with mean f i and variance S for the i-th time lag, i.e., t i ; N acf is the total number of ACF data points. Optimization of the logarithm of L acf with respect to provides the estimation of with the corresponding standard deviation. Note that the least-squares fitting method also works on the same assumption and approaches to same results 38 . However, we observe that such fitting predominantly depends on the range of the data considered. In Tab. 1, we show the estimated parameters with one standard deviation by fitting the ACF measured from a simulated data of total sampling time T s /τ 1 = 2000 and sampling time step �t/τ S = 0.1 with the same parameters as described above as input parameters. We vary the fitting range from time lag 0 − 0.1 upto 0 − 30 and show that the estimated results are closest to the input parameters for fitting rang 0-0.5 to 0-1. Note that in this range the standard errors corresponding to all the parameters are also the smallest.
To better understand the discrepancies in the obtained parameter values for different fitting ranges, we plot the ACF with the corresponding best fit (where the standard deviations are least) and residue in Fig. 2a (from simulated data of T s /τ 1 = 2000 and �t/τ S = 0.1 ) and Fig. 2b (from simulated data of T s /τ 1 = 10000 and �t/τ S = 0.1 ). Vividly, the central assumptions of the method, i.e., Gaussian distribution around the trend Table 1. Parameter estimated from simulated time series of T s /τ 1 = 2000 and time step �t/τ S = 0.1 by fitting the corresponding ACF over various range of time lags. Clearly, the estimations are close to all the input parameters when the fitting range lies in between from 0-0.5 to 0-1 s. The input parameters in the simulation are the same as in Fig. 1.

Fitting range (s)
η * 0 ± 1σ (mPa s) η * 1 ± 1σ (mPa s) τ The best fitting range increases with the increase in the total sampling time. Therefore, to make this process consistent, further processing of the data, such as careful binning is required 41 . Note that binning may lose information hidden in short times. If the sampling time step is much smaller compared to the smallest time scale τ S of the system then the binning can improve the fitting significantly. The parameter values obtained from the best fitting of the ACF are very close to those calculated using the Bayesian method developed in this paper (as shown in Fig. 1b-d).

Experiment
As viscoelastic fluids, we use an aqueous solution of polyacrylamide (Polysciences, PAAM, M w = 18 × 10 6 gm mol −1 ) at different concentrations. Note that the concentrations of the PAAM strongly dictates the viscoelastic parameters of the solution 25,42 . A very small volume fraction ( ∼ 0.01 ) of spherical polystyrene particles (radius 1.95µ m) are added to the PAAM-water solution as probes. For sample making, we seal the suspension between two glass slides with a spacing of 100 µ m. We kept this sample cell in contact with a thermal bath kept at 25 • C. We trap one of the spherical probe particles in the sample fluid far away from any glass surface (to avoid wall effects) using an optical tweezers built by tightly focused a Gaussian infra red (IR, = 1064 nm) laser beam with a high numerical aperture ( N.A. = 1.3 ) oil immersion microscope objective. We use a CMOS camera to record videos for the thermal position fluctuation of the probe particle at a frequency of 500 Hz. The position of the particle are tracked from the videos using MATLAB by standard particle tracking algorithm. For optimization of the log-likelihood function, we use an inbuilt MATLAB optimization tool. In Fig. 3a, we show a schematic of the setup, and Fig. 3b-c represent a typical position time series of the probe particle and the corresponding probability distribution respectively.

Data and analysis
For all measurements in this work, we maintain the trapping laser power fixed which results in trap stiffness k = (1 − 1.6) µN/m. Note that, the PAAM concentration may slightly change the refractive index of the fluid, which can bring a small modification to the trap stiffness. The trap stiffness is kept relatively lower to avoid any potential coupling of the trap with the fluid 43 . To experimentally check the dependency of the Bayesian method on the data length and the sampling time step, we analyze data for sampling times ranging from 100 to 4000 s and sampling time steps t = 0.002 s and 0.02 s. We show representative data ( 0.06 % w/w PAAM-water) in Fig. 4a-d for the estimations of η 1 τ 1 , η 0 and k respectively. Clearly, after T s = 200 s, the estimations are saturating for both the values of t . However, the precision and accuracy of the measurements are higher for t = 0.002 s. The standard error of the estimations, as expected, decreases with increasing T s . Interestingly, the best fittings of the auto-correlation functions evaluate approximately the same values as obtained from the Bayesian method, and exhibits similar behavior. As expected, the experimental data also shows the same strong dependence of the ACF fitting method on the range of data used for the fitting, as observed for numerically obtained ACFs. This verifies our theoretical understandings about the fitting method as described in the theory section. In Fig. 5d, we show the residue www.nature.com/scientificreports/ (inset) of the best fitting which clarifies that the assumption of the Gaussian distribution around the trend and the constant variance is only valid over the initial range of time lags. In order to verify the estimated results obtained from the Bayesian method, we perform the experiments with various (0.02-0.08% w/w) PAAM-water solutions and compare with the parameters obtained from the bulk rheology using a commercial rheometer (DHR-2-TA Instruments). We fit the frequency dependent viscosity measured using the rheometer with the corresponding theoretical expression for Jeffrey's fluid 24 to estimate the parameters. We observe (Fig. 5a-c) a very good agreement for all the three methods (Bayesian, best fitting of the ACF and rheometer measurements). The η 0 remains close to the viscosity of water at 25 • C and approximately independent of the concentrations, however, τ 1 and η 1 increases with the increase of the PAAM concentration. The similar trend is also observed previously in several works 25,42 . The standard errors are calculated over several independent measurements.

Conclusion
In conclusion, we present a Bayesian inference based technique to measure the rheological parameters of a Jeffrey's fluid. We show that the estimations from the proposed technique are more reliable compared to that obtained from the auto-correlation function. The parameter estimations from the ACF depend very sensitively on the range over which the data is fitted. We define the best fitting when the standard deviations corresponding to the estimations are minimum. At the best fitting, the estimations from the fitting method and the proposed method are approximately same. Further, to verify our method, we compare our inferred parameter values for various viscoelastic fluids prepared by mixing PAAM into water in different concentrations, with that measured using a commercial rheometer and obtain a strong agreement. We also observe that for a reliable measurement, the data should have sampling time step comparable to the shortest time scale of the system and total sampling time larger than ∼ 200 times the time constant of the fluid. The established method is comparatively fast (as it does not need further processing of the data), advanced, more reliable, and less affected by systematic noises. It uses less inputs from the user than any other least-squares fitting methods as it efficiently utilizes the information hidden in the data. Even though the Jeffrey's fluid model is a particular simple description of viscoelastic fluids, there exist many examples where such approach has been experimentally confirmed [8][9][10][11][12][13][14][15] . It should be mentioned,  www.nature.com/scientificreports/ however, that the concept of Bayesian inference can be also applied to more complex models involving more than a single stress relaxation time (even though the accuracy of the method generally becomes worse progressively with the introduction of additional parameters). The standard Bayesian model comparison techniques can be used to choose an appropriate model for a given data to avoid unnecessary generalization. Also, it can be extended for free Brownian probes and active Brownian particles. License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.