Microrheometer and viscoelastic properties of cells

The characterization of the mechanical properties of the cell is of great importance for estimating its biological state and physiology. Among the different measurement techniques, through the microrheometer it is possible to characterize the viscoelastic properties.

In this report we consider a measurement experiment of a paramagnetic bead microrheometer and the fitting of experimental data to analyze the behavior of a mechanical cell model with lumped parameters.

This model is described by the parallel of a Maxwell body with an elastic stiffness to which a damper in series is added to reduce the discrepancy with the experimental data. The viscoelastic constants estimated by the lumped parameter model are indicative of the properties of the entire cell and not of the individual constituents such as membrane, cytoskeleton or cytoplasm, for which a different model would be required.

The following report is an extract from a project carried out for the course of Modeling and Simulation of Physiological Systems, held for Medical Engineering at the University of Rome Tor Vergata.

Introduction

To mechanically characterize the cell and its viscoelastic properties it is necessary to have a descriptive model of its behavior. It is therefore considered a lumped parameter model consisting of a damper in series to a branch given by the parallel of a Maxwell body and an elastic stiffness. We therefore consider material parameters known from the literature through which the analytical and numerical solution to different types of input is analyzed.

Background

Cellular properties and physiology are the main reason for investigation and research. Among these mechanical properties are useful biomarkers [1].

The deformability and viscoeleastic properties are useful indicators of cytoskeletal changes and cellular state. Examples are the increase in stiffness in red blood cells affected by spherocytosis [2, 3] and the increase in deformability in cancer cells [4]. These mechanical properties can be measured with different techniques such as atomic force microscope, optical traps, magnetic ball microrheometer and aspiration with micropipettes [5].

This report focuses on a magnetic bead microrheometer and the cell model to estimate its viscoelastic properties.

Magnetic bead microrheometer

The microrheometer is a device used in micro rheology to characterize the rheological properties of a medium such as viscosity or flow traces. These devices are measured and are divided into active and passive depending on whether they exploit the thermal energy or external forces applied to move the particle.

In this report a magnetic bead microrheometer is analyzed. These 4.5 𝜇𝑚 diameter paramagnetic beads are linked to the cell cytoskeleton through the connection between the membrane integrins and the fibronectin which the beads are coated with. Membrane integrin is directly connected to the cytoscehelter so the measurement of the displacement of the ball provides information on cellular mechanics.

Magnetization

The presence of a controlled magnetic field applies a force to the marble and by measuring its displacement it is possible to characterize some viscoelastic properties [6] by fixing a viscoelastic model of the cell. The observed behavior is represented in fig. 1b and it is possible to distinguish three main phases. As the force is applied there is a first movement of the ball (phase I) identified as a fast elastic response. The first phase is followed by the relaxation regime (phase II) and finally the flow regime, where a progressive slip occurs.

Considering the following description of the model, it is possible to experimentally estimate three parameters such as the amplitude of the first phase 1 / 𝑘1 + 𝑘0, the time constant of the exponential trend of phase II equal to 𝜏 = 𝛾1 (𝑘0 + 𝑘1) / 𝑘0 𝑘1 and the parameter 𝛾0 identifying the amplitude of phase III. These values, estimated by Bausch et al. [6], will be used to implement a numerical solution of the model within the Matlab environment.

Mechanical model

These experimental curves are analyzed on the basis of a mechanical model described by a series of a damper (hereinafter referred to as Dashpot) and a Zener body, also known as the Kelvin body, described by a parallel of a branch with an elastic stiffness and a second branch formed by the series of a damped and a stiffness. This model is represented in fig. 1a. The mechanical description follows from the fact that the resulting displacement will be the sum of the dashpot 𝛾0 (D) and the Zener body (Z).

\begin{equation}
x(t)=x_{D}(t)+x_{z}(t)
\end{equation}

And therefore deriving it still holds:

\begin{equation}
\dot{x}(t)=\dot{x}_{D}(t)+\dot{x}_{z}(t)
\end{equation}

Clearly the force applied will be the same on both.

For (D) the relation of the damper is known from which we can derive the speed as:

\begin{equation}
\dot{x}_{D}=\frac{F}{\gamma_{0}}
\end{equation}

For the Zener body the force is distributed on the elastic branch and on the branch where there is a Maxwell body (spring and damper).

\begin{equation}
F(t)=F_{\text {Maxwell }}(t)+F_{\text {Ko }}(t)
\end{equation}

Where the force applied on Maxwell’s body is precisely the difference between the total applied force and that of the elastic branch 𝑘0𝑥𝑍 (𝑡). That is, for the Maxwell body, deriving the displacement, we obtain:

\begin{equation}
\dot{x}_{Z}(t)=\frac{\dot{F}}{k_{1}}+\frac{F}{\eta_{0}}
\end{equation}

And therefore the following applies to the Zener body as a whole:

\begin{equation}
\dot{x}_{Z}(t)=\frac{\dot{F}(t)-k_{0} \dot{x}(t)}{k_{1}}+\frac{F(t)-k_{0} x(t)}{\eta_{0}}
\end{equation}

In conclusion, the displacement will be given by the sum of the two displacements that can be obtained by solving the two differential equations for a given input and initial condition.

Results

This model is then analyzed within the Matlab numerical environment by comparing the numerical and analytical solution for different types of inputs.

Step input

For a step input, of amplitude ̄ 𝐹 we will obtain the differential equation for the dashpot:

\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{D}=\frac{\bar{F}}{\gamma_{o}} \quad t>0 \\
x_{D}(0)=0
\end{array}\right.
\end{equation}

Linear analytical solution over time:

\begin{equation}
x_{D}(t)=\frac{\bar{F}}{\gamma_{0}} t
\end{equation}

This linear response is representative of the linear phase observed in the experimental measurements. For the Zener body the following applies:

\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{z}(t)\left(1+\frac{k_{o}}{k_{1}}\right)=-\frac{k_{o}}{\gamma_{1}} x_{z}+\frac{\bar{F}}{\gamma_{1}} \quad t>0 \\
x_{z(0)}=\frac{\bar{F}}{k_{1}+k_{o}}
\end{array}\right.
\end{equation}

Where the initial condition derives from the fact to consider the displacement on the damper null at the initial instant and then the overall displacement is linked only to the elastic stiffnesses.

function dydt=zener_displacement(t,y,parameters,tspan,flag)
k_0=parameters(1); 
k_1=parameters(2);
gamma_0=parameters(3); 
gamma_1=parameters(4);
F_bar=parameters(5); 
switch flag
case ’square’
[F, dF]=force(t,tspan);
dydt=((F-k_0*y)/gamma_1 + (dF/k_1))*(1+k_0/k_1); 
case ’step’
dydt=(-(k_0/gamma_1)*y+(F_bar/gamma_1))/(1+(k_0/k_1)); 
case ’harmonic’
omega=parameters(7); 
F=F_bar*(1+sin(omega*t));
dF=F_bar*omega*cos(omega*t); 
dydt=((F-k_0*y)/gamma_1 + (dF/k_1))*(1+k_0/k_1);
end 
end

function dydt=dashpot_displacement(t,y,parameters,tspan,flag) gamma_0=parameters(3);
gamma_1=parameters(4); 
F_bar=parameters(5);
switch flag 
case ’square’
[F, dF]=force(t,tspan); 
dydt=F/gamma_0;
case ’step’ 
dydt=F_bar/gamma_0;
case ’harmonic’ 
omega=parameters(7);
F=F_bar*(1+sin(omega.*t)); 
dF=F_bar*omega*cos(omega*t);
dydt=F/gamma_0;
end
end

Considering the relaxation time:

\begin{equation}
\tau=\gamma_{1} \frac{k_{1}+k_{1}}{k_{0} k_{1}}
\end{equation}

We can rewrite the problem as:

\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{z}(t)+\frac{1}{\tau} x_{Z}(t)=\frac{\bar{F}}{k_{0} \tau} \quad t>0 \\
x_{z(0)}=\frac{\bar{F}}{k_{1}+k_{o}}
\end{array}\right.
\end{equation}

That is, this Cauchy problem presents a solution such as:

\begin{equation}
x_{z}(t)=e-{ }^{A(t)}\left[x_{0}+\int_{t_{0}}^{t} g(x) e^{A(t)}\right]
\end{equation}

Where 𝐴 (𝑡) is the term linked to the solution of the homogeneous equation, that is the integral of the exponential with the exponent the coefficient that multiplies 𝑥𝑍 (𝑡), while the term 𝑔 (𝑥) is linked to the non-homogeneous equation. Then we get:

\begin{equation}
x_{z}(t)=e-\frac{t}{\tau}\left[\frac{F}{k_{0}+k_{1}}+\frac{F}{k_{0}} e^{\frac{t}{\tau}}-\frac{F}{k_{0}}\right]
\end{equation}

Which can be rewritten as:

\begin{equation}
x_{Z}(t)=\frac{\bar{F}}{k_{0}}\left(1-\frac{k_{1}}{k_{0}+k_{1}} e^{-\frac{t}{\tau}}\right)
\end{equation}

That is an exponential that will progressively tend to ̄ 𝐹 / 𝑘0.

Furthermore, the displacement is depicted in fig. 3a and fig. 9. This solution is graphed in fig. 3a and compared with the numerical solution obtained in Matlab.

Numerical solution

The numerical solution is obtained in Matlab through 𝚘𝚍𝚎𝟷𝟻𝚜, a proprietary solver for stiff differential equations.

The main routine (main.m) allows you to select the type of force input using the variable 𝚏𝚕𝚊𝚐, choosing between the step input, a square wave and a sine wave, described below. Subsequently the routine sets an appropriate time interval and calls 𝚘𝚍𝚎𝟷𝟻𝚜 () passing it the functions where the displacements are defined (fig. 5) and then their sum is plotted.

Square wave input

A square wave input is then tested as used in [6], with a period of 5s, a duty cycle of 50% and an amplitude of 2000 pN, tested over an interval of 20s. This input is implemented through the force.m routine, described in fig. 8 having the foresight to soften the step in order to be able to make the numerical derivative even of the discontinuities and making sure that 𝚘𝚍𝚎𝟷𝟻𝚜 can request the value at any instant, and not only on the sample where it is defined, through an interpolation.

This signal is present in fig. 7. The results are shown in fig. 3b and from the decomposition of the two signals (fig. 9) it is possible to observe how the dashpot has a piecewise linear component to which the oscillatory behavior of the Zener body is added.

Sinusoidal input

The oscillatory effect is clearly very present if we stimulate the model with a sinusoidal input of the type:

\begin{equation}
F(t)=\bar{F}(1+\sin (\omega t))
\end{equation}

The analytical solution has both a harmonic and an exponential component:

\begin{equation}
\begin{aligned}
x(t)=\frac{F_{0}}{k_{0}} & {\left[1-\frac{k_{1}}{k_{0}+k_{1}} e^{-t / \tau}+\right.} \\
&+\frac{\tau \omega\left(1-\frac{\gamma_{1}}{k_{1} \tau}\right)\left(e^{-t / \tau}-\cos \omega t\right)}{1+(\tau \omega)^{2}}+\\
&\left.+\frac{\left.\left(1+\frac{\gamma_{1} \tau \omega^{2}}{k_{1}}\right) \sin \omega t\right]}{1+(\tau \omega)^{2}}\right]
\end{aligned}
\end{equation}

This solution is compared with the numerical solution in fig. 6. From the two separate displacements (fig. 9) it is possible to see how both have a harmonic component and while the displacement of the Zener body is limited, the dashpot presents a linear growth enveloped by a sinusoid.

function [F, dF]=force(t_curr,t)
 T_cyc =5;
 N_cyc =4;
 F=2000;
 freq=1/T_cyc;
 squareWave=F/2*(square(2*pi*freq*t,50)+1);
 squareWave(end)=0;
 squareWave_smooth=smoothdata(squareWave ,’gaussian’,(size(t)/N_cyc)*0.05);
 derivate=[0 diff(squareWave_smooth)];
 F=interp1(t,squareWave_smooth ,t_curr);
 dF=interp1(t,derivate ,t_curr);
end

Conclusions

The study and classification of the mechanical behavior of these models makes it possible to analyze their coefficients and to be able to use experimental data to derive the cellular mechanical properties. Using a numerical solver allows you to analyze the mechanical model even for complex and highly variable inputs over time.

The viscoelastic constants estimated by the lumped parameter model are indicative of the properties of the entire cell, but these properties are determined in part by the properties of the cell membrane and in part by the properties of the cytoplasm and the structure of the cytoskeleton. A more detailed model is required for more information on the individual cellular components.

Code availability

The code is available at the project’s GitHub repository.

References

• [1] Dino Di Carlo. “A Mechanical Biomarker of Cell State in Medicine”. en. In: SLAS Technology 17.1 (Feb. 2012), pp. 32–42. ISSN: 24726303. DOI: 10 . 1177 / 2211068211431630.
• [2] Gabriel Y.H. Lee and Chwee T. Lim. “Biomechanics approaches to studying human diseases”. en. In: Trends in Biotechnology 25.3 (Mar. 2007), pp. 111–118. ISSN: 01677799.
• [3] S. Suresh et al. “Connections between single-cell biomechanics and human disease states: gastrointestinal cancer and malaria”. en. In: Acta Biomaterialia 1.1 (Jan. 2005), pp. 15–30. ISSN: 17427061.
• [4] Claudia Tanja Mierke. “Viscoelasticity Acts as a Marker for Tumor Extracellular Matrix Characteristics”. In: Frontiers in Cell and Developmental Biology 9 (Dec. 2021), p. 785138.
• [5] C Ross Ethier and Craig A Simmons. Introductory Biomechanics: From Cells to Organisms. en.
• [6] Andreas R. Bausch et al. “Local Measurements of Viscoelastic Parameters of Adherent Cell Surfaces by Magnetic Bead Microrheometry”. en. In: Biophysical Journal 75.4 (Oct. 1998), pp. 2038–2049. ISSN: 00063495.