Recursive Least Squares Filtering Algorithms for On-Line Viscoelastic Characterization of Biosamples

The mechanical characterization of biological samples is a fundamental issue in biology and related fields, such as tissue and cell mechanics, regenerative medicine and diagnosis of diseases. In this paper, a novel approach for the identification of the stiffness and damping coefficients of biosamples is introduced. According to the proposed method, a MEMS-based microgripper in operational condition is used as a measurement tool. The mechanical model describing the dynamics of the gripper-sample system considers the pseudo-rigid body model for the microgripper, and the Kelvin–Voigt constitutive law of viscoelasticity for the sample. Then, two algorithms based on recursive least square (RLS) methods are implemented for the estimation of the mechanical coefficients, that are the forgetting factor based RLS and the normalised gradient based RLS algorithms. Numerical simulations are performed to verify the effectiveness of the proposed approach. Results confirm the feasibility of the method that enables the ability to perform simultaneously two tasks: sample manipulation and parameters identification.


Introduction
The mechanical characterization of biomaterials represents a crucial procedure in the fields of tissue engineering, tissue and cell mechanics, disease diagnosis and minimally invasive surgery (MIS) [1][2][3].
For instance, modeling the mechanical response of brain tissue can be useful in the understanding of traumatic brain injury [4] or blast-induced neurotrauma [5] mechanisms.In tissue engineering, scaffolds should achieve appropriate mechanical properties in order to ensure tissue growth.Therefore, before implantation, it is necessary to verify if these properties meet the specified requirements [6,7].Many investigations have also been focused on the characterization and on the constitutive modeling of skin [8,9], cornea [10], vocal fold [11], skeletal muscle [12], and blood vessels [1] tissues.Regarding the circulatory system, it is worth mentioning that vascular stiffness is recognized as an important factor in diseases such as pulmonary arterial hypertension, kidney disease, and atherosclerosis [13].Generally, the variation of tissue stiffness is a hallmark of several disease states, including fibrosis and some types of cancers [14].Viscosity can also be considered as a biomarker for the metastatic potential of cancer cells [15,16], and the analysis of both elastic and viscous properties could be more effective in detecting and identifying specific diseases [17].
The characterization of the mechanical properties of biomaterials is important not only for in vitro analysis.For example, in MIS operations, tactile information is not available to the surgeon.Therefore, the capability of on-line detecting and identifying tissues could play a fundamental role in the advancement of MIS procedures [3].
The mechanical properties of tissues have been investigated by using various methods, such as indentation and aspiration, and measuring deformations by means of ultrasound or magnetic resonance imaging techniques [18].Surgical instruments have also been improved with sensing capabilities [19].For example, a micro-tactile MEMS-based sensor to be integrated within MIS graspers was presented in Reference [3], whereas a sensor capable of determining the stiffness of biological tissues was modeled and tested in Reference [20].
At the cell level, experimental techniques for investigating mechanical properties include force-application techniques (e.g., micropipette aspiration, atomic force microscope probing, optical trapping), and force-sensing techniques (e.g., traction force microscopy, wrinkling membranes, micropost arrays).The technique selection depends on size of the biological sample, feature to be acquired, detrimental effects on the sample, spatial and force resolution, and accuracy [21,22].
Given the large number of experimental approaches, many mechanical models have been proposed in literature, based on the micro/nanostructural approach and on the continuum approach.The former model targets the sub-cell level, whereas the latter one treats the biological sample as a continuum material.Considering the continuum approach, many constitutive models have been proposed in literature, such as biphasic model, liquid drop models, solid viscoelastic models (Kelvin-Voigt, Maxwell, Zener) [23,24].Once the mathematical model is available, a parameter estimation problem can be formulated and solved by means of different methods, such as least-squares [25], Kalman filtering [26], or inverse finite element [18,27] methods.
In the last decade, microgrippers became essential tools in the manipulation at the microscale [28][29][30][31][32][33], and can be adopted to develop novel techniques for the viscoelastic characterization of soft materials.In their previous investigation [34], the Authors proposed a method to evaluate the mechanical properties of a biomaterial sample based on the use of a MEMS microgripper.The mechanical model was developed considering the gripper-sample kinematics and the Kelvin-Voigt constitutive law of viscoelasticity, and a PID controller to determine stiffness and viscous parameters of the sample.
In this work, a parameter estimation problem is formulated starting from the model presented in Reference [34].Recursive least squares filtering algorithms are implemented to characterize the mechanical properties of soft biomaterials, enabling the on-line identification of stiffness and viscosity coefficients.With respect to the other experimental techniques briefly mentioned above, that generally are implemented only for the mechanical testing procedure, the proposed approach enables the ability to perform the mechanical characterization while the sample is manipulated for a different purpose in another operation.Furthermore, this method does not require additional elements other than the measurement tool and the specimen as, for example, beads used with optical tweezers.Also, risks related to laser-induced damage and to large deformations of the sample (involved in optical traps and micropipette aspiration, respectively) are avoided.
The paper is organised as follows.The microdevice is introduced in Section 2, whereas the mathematical model is discussed in Section 3. The adopted identification algorithms are described in Section 4, and simulations and numerical results are reported and discussed in Section 5.

The Experimental Device
The microsystem considered in this paper consists of two silicon arms actuated by rotary comb drives, as shown in Figure 1.The electrostatic actuators exert the input torques that, through the conjugate surface flexure hinges (CSFH) deflections [35,36], move the arms from the initial position to perform the gripping task [37].In the neutral configuration, the distance between the jaws is equal to 150 µm.The fabrication method relies on the application of deep reactive-ion etching on silicon-on-insulator wafers [38].In their previous investigations [34,39], the authors estimated the mechanical characteristic of the sample considering input signals of suitable waveform.According to the proposed method, the sample is gripped with the left arm until the right arm reaches a predefined rotation, that serves as a reference signal for the implemented feedback control scheme.Therefore, the stiffness coefficient can be computed at steady state conditions [39].In order to estimate also the sample viscosity, a small-amplitude sinusoidal signal was added to the left arm, and the viscous coefficient was obtained as a function of the input torque frequency [34].
In this paper, an estimation algorithm for the elastic and viscosity parameters of the gripped sample under generic control torques is presented.As demonstrated in the next Sections, the algorithm take advantage of the mathematical model structure that, despite its general non linearity, is linear with respect to the unknown parameters.

The Mathematical Model
In operative condition, the microdevice is gripping a sample between its jaws, as illustrated in Figure 2.More specifically, Figure 2a shows a schematic drawing of the compliant structure of the gripper.The contact points between jaws and sample are labelled as B and C, whereas the points A and D represent the centers of rotation of the left and right arms, respectively.Starting from the compliant mechanism, a pseudo-rigid body model (PRBM) can be obtained by substituting the constant-curvature CSFHs with revolute joints [35,[40][41][42].Figure 2b shows the PRBM of the microgripper: the quadrilateral ABCD is a closed chain composed by the links AD (frame), AB and CD (arms), and BC (sample).All the links have a fixed length with the exception of BC, whose compression is at the basis of the measurement procedure.The proposed model takes into account the nonlinearity related to the kinematics of the gripper, but it assumes a linear behavior for stiffness coefficients k 2 and k 4 .Moreover, the identification of the cell parameters inherently assumes that its dynamic behavior can be described by the Kelvin-Voight constitutive model of viscoelasticity.
With reference to Figure 2b, the following parameters can be introduced: • l is the common length of the two links representing the arms, i.e., the distances AB and CD; • d is the distance between the hinges (AD), i.e., the frame length; • k 2 , k 4 and k are the torsional stiffness of the two arms hinges and the stiffness of the sample, respectively; For the sake of simplicity in the model representation, all the considered variables are referred to the neutral layout: the gripper, in symmetrical configuration, is in contact with the sample but no deformation occurs.By following the notation introduced in Reference [34], and with reference to Figure 3a, the angles θ2 , θ3 , and θ4 represent the reference orientations of the links AB, BC and DC, respectively, whereas û represents the reference length of the link BC.The notation • is used for the parameters in the deformed configuration, as shown in Figure 3b.Therefore, the angles θ i = θi − θi represent the relative angular displacements of the two links from their neutral configuration, with i = 2 for the left link and i = 4 for the right one.The orientation of BC follows the same notation, with i = 3 and θ3 = 0, whereas the deformation is equal to u = ũ − û.The values of the variables in the neutral configuration are given in Table 1.Assuming the inertia of the sample to be negligible, the dynamical model of the links can be written as where the subscripts 2 and 4 refer to the left and right joints, respectively.For the device here considered, the values of the parameters appearing in ( 1) and ( 2) are reported in Table 2.
The system has two degrees of freedom and it is fully described by Euqations (1) and ( 2).The variables θ 3 and u, and their time derivatives, can be computed as functions of the state variables θ 2 , θ 4 , θ2 and θ4 [34]:

Mechanical Characteristics Estimation
As briefly described in the Introduction, an approach to determine both the elastic and viscous characteristics of a biosample by using a microgripper was proposed in Reference [34], where particular waveforms were required for the input torques.However, this method is not suitable to perform simultaneous manipulation and measurement tasks.
To overcome this limitation, a measurement scheme adaptable to various operative conditions should be considered.Therefore, an on-line dynamical estimator could be implemented as a more efficient strategy for the mechanical characterization of the sample.Recursive least square methods can be successfully applied when the unknown parameters are the linear coefficients of the system dynamic equations.In such cases, quite usual in literature [43][44][45][46], the model can be rearranged as a linear time-varying system with respect to the parameters to be estimated, whereas the other terms, supposed measurable, are functions either of the state or of the output variables.
The model dynamics has to be rearranged as where m is the number of degrees of freedom of the system (m = 2 in the considered case), ω i (t) is the vector of unknown parameters, and y i (t), M i (t) are known quantities depending on the system dynamics.It is worth noting that all the terms in Equation ( 6) are time dependent.Equations ( 1) and ( 2) can be rewritten, according to the notation in Equation ( 6), in linear form with respect to the parameters as where the sign in Equation ( 8) is minus or plus if i = 2 or if i = 4, respectively.
Referring to Equations ( 7)-( 9), the general expressions to be defined for a generic recursive last squares (RLS) filtering algorithm are [47]: where ωi (t) and ŷi (t) are the current estimation values of ω i (t) and y i (t), i (t) is the current prediction error, K i (t) is the gain determining how much the prediction error affects the update in the parameters estimation, and φ i (t) represents the gradient of the predicted model output with respect to ω i (t).
Two recursive last squares (RLS) filtering algorithms are applied considering different ranges of values for the parameters to be estimated.The algorithms are then compared to each other to understand how much the different viscous and elastic characteristics of the dynamical system affects the convergence of the estimation algorithms.With respect to the formulation in Equation ( 10), the two approaches differ for the choice of Q i (t).Moreover, it is worth noting that the symmetric structure of the dynamics represented by Equations ( 1) and (2) implies that, under a full state measurement, the results obtained choosing i = 2 are equal to the results corresponding to the case i = 4. Therefore, firstly the case i = 2 and consecutively the case i = 4 are addressed.Once equivalence and effectiveness of both cases has been proved, the choice in real applications should be driven by the simplicity of the measurement.For example, the operative technique of torque action on the first jaw only implies that for i = 4 no torque measurement is required, strongly simplifying the implementation.

Forgetting Factor Based RLS
The first method is based on a forgetting factor based RLS algorithm.In Equation (10), the following choice is performed where and It's assumed that the residual i (t) (the difference between the estimated and the measured value of y i (t)) is affected by a white noise with covariance equal to 1.According to previous equations, ωi (t) is computed in order to minimize the sum of residuals squares ωi (t) = arg min In Equations ( 11)-( 14), λ ∈ R is the forgetting factor introduced in order to consider differently the time sequence of the errors i (t), according to an exponentially decreasing weight if λ ∈ (0, 1).This choice is effective in case of time varying parameters while, when dealing with constant parameters, the choice λ = 1 is usually adopted.
The algorithm in Equation ( 10), with positions ( 11)-( 13), has been applied to the case i = 2.In all the simulations, the initial values of the parameters have been chosen far from the real values.The initial covariance, proportional to P 2 , has been fixed taking into account that the covariance matrix has to be chosen according to a priori knowledge of the parameters at t = 0: very high values of the covariance matrix elements correspond to completely unknown parameters.
Remark 1.Note that the forgetting factor method is a particular simplified case of the Kalman filter.

Normalised Gradient Based RLS
The second approach consists of a normalized gradient algorithm.This technique is based on the choice of Q i (t) with a simpler form with respect to (11): where γ is the adaptation gain scaled by the gradient φ i (t).The bias term β is added to the square norm of the gradient vector in the denominator of the expression, in order to prevent critical situations in the case φ i (t) is close to zero.The algorithm (10) with Q i (t) defined by (15) requires only the initialization for the values of the parameters to be estimated.Since the presence of noise is not explicitly considered in this formulation, the drawback is a smaller rate of convergence and a larger sensitivity to the presence of noises.
This technique has been simulated making reference to the dynamics represented by Equation (2), that is with i = 4.

Simulations
Numerical simulations, using Matlab R (MathWorks, Inc., Natick, MA, USA) and Simulink R (MathWorks, Inc., Natick, MA, USA) tools, were performed in order to show effectiveness, benefits and differences of the proposed estimation methods.Three different numerical cases were analysed, considering predominant elastic or dissipative behaviours.
The first case corresponds to a realistic condition with elastic and damping coefficient much greater than the ones of the mechanical structure, with c = 8.4 × 10 −6 Nms/rad and k = 2.5 × 10 −3 Nm/rad, where the elastic coefficient greater than the damping one.
In the second case, a damping coefficient greater than the elastic one was considered in order to check, by comparison, the dependency of the algorithm convergence from the two different mechanical characteristics.The order of magnitude for the two coefficients have been exchanged, setting c = 8.4 × 10 −3 Nms/rad and k = 2.5 × 10 −6 Nm/rad.
The initial parameters values, for all the simulations and for both the estimation algorithms have been chosen as c(0) = 10 −9 Nms/rad and k(0) = 10 −7 Nm/rad.
For the forgetting factor RLS estimator introduced in Subsection 4.1, the 2 × 2 square covariance matrix is diagonal, with both the diagonal elements equal to 10 20 , while the forgetting factor λ is fixed to λ = 0.99.
For the normalized gradient estimator, presented in Subsection 4.2, the adaptation gain γ has been set as γ = 0.9 and the normalization bias has been chosen as β = 2.2 × 10 −16 .
Simulation results obtained for the first case (c = 8.4 × 10 −6 Nms/rad and k = 2.5 × 10 −3 Nm/rad) are depicted in Figure 4a for the elastic coefficient k and in Figure 4b for the damping coefficient c.The solid line shows the estimation evolution with the first algorithm, whereas the dashed one reports the results of the second algorithm.The dotted line corresponds to the true values of the parameters, plotted as a reference.As expected, both the algorithms converge in a very short time, but the first is faster than the second.For the second case (c = 8.4 × 10 −3 Nms/rad and k = 2.5 × 10 −6 Nm/rad), the simulation results are depicted in Figure 5a for the damping coefficient c, and in Figure 5b for the elastic one k.As in the previous case, the solid line refers to the estimation evolution with the first algorithm, the dashed one refers to the second algorithm, and the dotted one is the true reference value.The difference in the convergence rate for the two approaches is confirmed.The results obtained by simulation of the third case (c = 8.4 × 10 −11 Nms/rad and k = 2.5 × 10 −5 Nm/rad) are reported in Figure 6a for the elastic coefficient k, and in Figure 6b for the damping one c.The difference in the convergence rate for the two approaches is confirmed, and uniformity in the estimation of the damping coefficient c with the two algorithms is the same as in the first case.Numerical results support the effectiveness of the proposed method, tested considering different orders of magnitude for the parameters to be estimated.Clearly, in the present formulation, the state is assumed to be measurable, with possible additive Gaussian noise to model uncertainties and measurement noise.

Conclusions
In this paper, the possibility of using a microgripper for the identification of the mechanical properties of biomaterials was investigated.To describe the gripper-sample system dynamics, the pseudo-rigid body model and the Kelvin-Voigt constitutive law of viscoelasticity were considered for the microgripper and for the sample, respectively.The normalised gradient-based RLS and the forgetting factor based RLS algorithms were implemented to solve the parameters estimation problem.Three cases were analysed, assuming different ranges of values for the coefficients to be estimated.The simulation results confirmed the feasibility of the method: both the algorithms converged in a very short time, but the normalised gradient-based RLS algorithm resulted to be faster than the other one.Therefore, the proposed approach could enable the ability to simultaneously perform the manipulation of the biosample and the identification of its mechanical characteristics, i.e., stiffness and damping coefficients.However, the application of the proposed method to a variety of biological samples is limited by the assumption of the Kelvin-Voigt constitutive law of viscoelasticity.From the methodological point of view, further investigations will focus on the implementation of more complex constitutive models and on the design of high-performance algorithms, always attaining robustness with respect to noise and parameter uncertainties.

Figure 1 .
Figure 1.Optical microscope images of the device.The whole microgripper (a) and two overlapping frames showing the left arm in neutral (0 V) and actuated (28 V) configurations (b).

Figure 3 .
Figure 3. Nomenclature and model parameters in neutral (a) and general (b) configurations.

Table 2 .
Numerical values of the parameters.