This paper has been accepted for publication in IEEE Transactions on Biomedical Engineering.

This is the author’s version of an article that has, or will be, published in this journal or conference. Changes were, or will be, made to this version by the publisher prior to publication.

DOI:	10.1109/TBME.2022.3179655
IEEE Xplore:	https://ieeexplore.ieee.org/document/9786847

Please cite this paper as:

A. Bianco, R. Zonis, A.-M. Lauzon, J. R. Forbes and G. Ijpma, “System Identification and Two-Degree-of-Freedom Control of Nonlinear, Viscoelastic Tissues,” in IEEE Transactions on Biomedical Engineering, doi: 10.1109/TBME.2022.3179655.

©2022 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

System Identification and Two-Degree-of-Freedom Control of Nonlinear, Viscoelastic Tissues

Amanda Bianco Raphael Zonis Anne-Marie Lauzon James Richard Forbes and Gijs Ijpma Manuscript received December 7, 2021; revised May 21, 2022; accepted May 22, 2022. This work was supported by the Natural Sciences and Engineering Research Council of Canada and Aurora Scientific Inc. through the Collaborative Research and Development grant CRDPJ-519929-17 for the Lauzon Lab, and by the Fonds de recherche - Nature et technologies Master’s Research Scholarship awarded to Amanda Bianco. (Corresponding author: A-M. Lauzon)A. Bianco was with the Department of Biomedical Engineering and Meakins-Christie Laboratories, McGill University, Montreal, QC, Canada. She is now with the Faculty of Medicine and Health Sciences, McGill University, Montreal, QC, Canada.R. Zonis was with the Department of Mechanical Engineering, McGill University, Montreal, QC, Canada. He is now with the Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA, USA.A-M. Lauzon is with the Department of Medicine, Department of Biomedical Engineering and Meakins-Christie Laboratories, McGill University, Montreal, QC, Canada (e-mail: [email protected]).J. R. Forbes is with the Department of Mechanical Engineering, McGill University, Montreal, QC, Canada.G. Ijpma is with the Meakins-Christie Laboratories, McGill University, Montreal, QC, Canada.This article has supplementary downloadable material available at https://doi.org/10.1109/TBME.2022.3179655, provided by the authors.Digital Object Identifier 10.1109/TBME.2022.3179655

Abstract

Objective: This paper presents a force control scheme for brief isotonic holds in an isometrically contracted muscle tissue, with minimal overshoot and settling time to measure its shortening velocity, a key parameter of muscle function. Methods: A two-degree-of-freedom control configuration, formed by a feedback controller and a feedforward controller, is explored. The feedback controller is a proportional-integral controller and the feedforward controller is designed using the inverse of a control-oriented model of muscle tissue. A generalized linear model and a nonlinear model of muscle tissue are explored using input-output data and system identification techniques. The force control scheme is tested on equine airway smooth muscle and its robustness confirmed with murine flexor digitorum brevis muscle. Results: Performance and repeatability of the force control scheme as well as the number of inputs and level of supervision required from the user were assessed with a series of experiments. The force control scheme was able to fulfill the stated control objectives in most cases, including the requirements for settling time and overshoot. Conclusion: The proposed control scheme is shown to enable automation of force control for characterizing muscle mechanics with minimal user input required. Significance: This paper leverages an inversion-based feedforward controller based on a nonlinear physiological model in a system identification context that is superior to classic linear system identification. The control scheme can be used as a steppingstone for generalized control of nonlinear, viscoelastic materials.

{IEEEkeywords}

Feedback, feedforward, force control, muscle.

1 Introduction

\IEEEPARstart

Muscle mechanics measurements provide key insights into the contractile and non-contractile properties of healthy and diseased muscles. Additionally, these measurements can be used to assess the response to drugs, biological mediators and/or environmental factors. While crude mechanics measurements can be performed in vivo, isolated muscle bundles need to be studied in vitro in order to obtain accurate control over contractile stimulation in addition to direct measurement and control of length and force. In vitro muscle mechanics measurements can be done at the cell or tissue level. At the cell level, techniques like traction force microscopy and collagen gel constructs can directly measure the forces that cultured muscle cells exert on a substrate [1, 2, 3]. However, studies show that individual muscle cells are prone to losing their contractile phenotype when studied in cell culture [4, 5, 6]. To better characterize in vivo muscle mechanics, freshly dissected, intact muscle tissue can be used directly as a higher fidelity model.

To measure the force generated by a muscle tissue, one end of the tissue is attached to a force transducer, allowing direct force readout. However, a number of key parameters of muscle contractility require direct control of the force. One such parameter is the muscle’s rate of shortening, otherwise known as the shortening velocity. The shortening velocity at no load, also known as the unloaded shortening velocity, is believed to be a direct correlate of the molecular level mechanics of the muscle’s contractile elements [7]. One way to measure the shortening velocity of muscle tissue in vitro is via quick-release experiments in which the tissue is suspended between a force transducer and a length actuator and subjected to either a chemical or an electrical stimulus to induce contraction [7]. During an isometric contraction of the muscle tissue, a sudden switch to isotonic force control is applied. The level of force to which the tissue is controlled will determine the rate at which the tissue shortens, i.e., the shortening velocity. As unloaded shortening velocity cannot be measured directly in an actively controlled tissue with no compressive strength, a series of different force level quick-release protocols, which are referred to as force clamps, are used to find the unloaded shortening velocity by extrapolating the hyperbolic curve that is fit to the data [8]. By the nature of the fit, the estimation of the unloaded shortening velocity is strongly dependent on force clamps closest to a purely unloaded state. In addition to acquiring shortening velocity measurements and an estimate of the unloaded shortening velocity, force clamps can also provide insight into other mechanical parameters such as elastic recoil. Force clamps were traditionally performed by attaching physical loads to tissue samples [9], but are now performed with simulated loads via digital force control [10, 11, 12, 13, 14].

The literature in the field of tissue mechanics does not provide much detail regarding the force control scheme used for force clamps since researchers often use commercial equipment with control schemes that are unpublished and require manual tuning [10, 11, 12, 13, 14]. The authors of [15] present the design of a feedback system involving proportional-integral (PI) control for isotonic studies of smooth muscle in vitro, but manual tuning is required. As a consequence of intra- and inter-tissue variability, coupled with the dynamic, nonlinear, viscoelastic nature of muscle tissue, it is difficult to design a static, linear feedback controller that can successfully perform force clamps in all test cases. This lack of versatility makes the performance of force clamps highly cumbersome, and precludes the throughput necessary to perform high volume studies, as well as the possibility of complete automation.

The contributions of this paper include leveraging a nonlinear physiological model in a system identification (ID) context that is superior to classic linear system ID and the implementation of an inversion-based feedforward controller to help automate force clamps for measuring muscle mechanics. While specifically applied to the determination of shortening velocity here, the approach can be generalized to any force control of nonlinear, viscoelastic materials. This paper is structured as follows: Section 2 describes the experimental setup and protocols needed to test the force control scheme, Section 3 details the force control scheme along with the control objectives and the system ID process, Section 4 presents the experimental results and Section 5 discusses the experimental results with respect to the control objectives.

2 Materials and Methods

2.1 Muscle Tissue Procurement and Preparation

The force control scheme was developed and tested in experiments with equine airway smooth muscle (ASM) and its robustness was confirmed with murine flexor digitorum brevis (FDB) muscle, which is a skeletal muscle. Equine trachea segments were obtained from a slaughterhouse, after which the ASM with epithelium and connective tissue was dissected out and cryostored according to [16]. On the day of the experiment, ASM tissues were dissected free from the epithelium and connective tissue in Ca²⁺-free Krebs-Henseleit solution (refer to Section 1.1 of the Supplemental Material for solution compositions). Aluminum foil clips were then attached to either end of each tissue, which was approximately 3 $\times$ 1 $\times$ 0.5 mm. The ASM tisues that were not immediately used were pinned to pieces of silicone at in situ length and kept in Dulbecco’s modified Eagle’s medium with 2 $\%$ heat-inactivated fetal bovine serum at 37^∘C in an incubator for no more than 3 days. Murine tissues were procured from waste tissue from mice used for unrelated research. All procedures with murine FDB muscle were approved by the McGill University Animal Care Committee and performed in accordance with the guidelines of the Canadian Council on Animal Care (protocol 8082 approved in June 2019). The FDB muscles were dissected from the hind paws according to [17, 18], removing as much of the surrounding tissue as possible. Aluminum foil clips were then attached to either end of each FDB muscle, which was approximately 5 $\times$ 1 $\times$ 1 mm. FDB muscles that were not immediately used were preserved in the same manner as ASM.

2.2 Tissue Bath Setup

2.2.1 Hardware

The tissue bath setup includes a force transducer (model 400A from Aurora Scientific Inc.) and length actuator (model 322C from Aurora Scientific Inc.) on either side of a tissue bath (Fig. 1). The resonant frequency of the force transducer and length actuator are approximately 140 Hz and 2 kHz, respectively, both of which are greater than the assumed bandwidth of muscle tissue, e.g., 100 Hz or less based on the frequency spectrum of the inputs explored in [19]. The voltage-to-force conversion factor of the force transducer was 5 mN/V and the voltage-to-length conversion factor of the length actuator was 0.3 mm/V. Fluid flow was controlled by two peristaltic pumps, one for relaxing solution and one for contracting solution, at 1 mL/min into a total bath volume of 2.4 mL, with overflow suction to a vacuum. Platinum electrodes, which were connected to an electrical stimulator (model 701C from Aurora Scientific Inc.), were placed in the tissue bath on either side of the muscle tissue for electric field stimulation (EFS).

Data acquisition, force control, peristaltic pump control and EFS control were achieved via a FPGA-controlled USB-7845 system from National Instruments.

Refer to caption — Figure 1: Tissue bath setup, where the muscle tissue is ASM.

2.2.2 Software

LabVIEW code, uploaded to the FPGA system, controlled the tissue bath’s hardware. MATLAB was primarily used for data analysis, but MATLAB scripts were executed within the LabVIEW environment of the host computer primarily for model estimation/validation and feedforward generation, all of which are described in Section 3.

2.3 Experimental Protocols

All experiments were conducted at room temperature $\pm$ 1 $\degree$ C. The tissue bath was filled with Ca²⁺-free Krebs-Henseleit solution and a muscle tissue was mounted between the force transducer and length controller (Fig. 1). Once the tissue was mounted, it was equilibrated with repeated contractions induced either by 10^-5 M methacholine solution for ASM or EFS for FDB muscle (refer to Section 1 of the Supplemental Material for more details on chemical and electrical stimulation parameters). After model estimation/validation and feedforward generation, force clamps at 6 levels, e.g., 5 $\%$ , 7 $\%$ , 10 $\%$ , 20 $\%$ , 40 $\%$ and 80 $\%$ of the steady state contractile force achieved directly prior to force clamp initiation, were applied in random order at 1-min intervals, all within one contraction for ASM and in separate contractions for FDB muscle on account of the short maximal duration of EFS. Note that 5 $\%$ and 80 $\%$ are the largest and smallest amplitude force clamp levels, respectively, since the reference contractile force was dropped by 95 $\%$ for the former and by 20 $\%$ for the latter. To test the repeatability of the force control scheme, force clamps were repeated within and across contractions of the same contractile force for ASM. As for FDB muscle, force clamps were repeated across contractions of the same as well as different contractile force. The repeatability and robustness of the force control scheme depends on that of the system ID process described in Section 3.3, specifically model estimation. Thus, it was important to assess whether the same model could be estimated each time this process was repeated and if the same estimate could then be used for each force clamp level within or across contractions of the same or different contractile force. This assessment was mainly done for the selected nonlinear model. Refer to Section 1 of the Supplemental Material for the details regarding the experimental protocols.

3 Force Control Scheme

3.1 Control Objectives

Different experiments have different requirements for settling time and overshoot based on the research needs. The maximum settling time is set here to 60 ms so that shortening velocity measurements can take place as early as 60 ms into the force clamp. The allowable overshoot is less than 3 $\%$ of the contractile force just prior to the force clamp, specifically for the larger amplitude force clamp levels, e.g., those of 10 $\%$ or less, to avoid the highly nonlinear tissue mechanics near slack length. The force control scheme also needs to be repeatable, robust and require little to no inputs and supervision from the user, which are especially important objectives for long-term automated tissue studies. It is important to note that shortening velocity is neither measured nor evaluated in this paper because there do not exist true shortening velocity values against which measurements can be compared. This is why the control objectives are focused on tracking the reference signal via the control parameters of settling time and overshoot, with the emphasis being on the larger amplitude force clamp levels for accurate extrapolation of the unloaded shortening velocity.

3.2 Control Configuration

3.2.1 Overview

The implemented two-degree-of-freedom (2DOF) control configuration [20] is shown in Fig. 2. It comprises the plant $P$ , feedback controller $C_{fb}$ and feedforward controller $C_{ff}$ . The difference between the reference signal $r$ and measured output $y$ is the error $e$

e=r-y,

(1)

which serves as the input to the feedback controller. The reference signal also serves as the input to the feedforward controller. The total control effort $u$ is given by

u=u_{ff}+u_{fb},

(2)

where $u_{ff}$ is the control effort of the feedforward controller and will be referred to as the feedforward signal and $u_{fb}$ is the control effort of the feedback controller. The feedforward and feedback controllers can be designed together or independently. In this work, for simplicity as well as to leverage system ID techniques and information provided by physiological models, they are designed independently.

3.2.2 Reference Signal

The reference signal consists of a 600 ms step to the desired force clamp level scaled by the reference contractile force. This reference force is determined by averaging 100 samples of force reading prior to initiating the force clamp. The reference signal is ramped down over 20 ms to help manage the hardware delay between the input and output of the plant as well as the resonance of the force transducer. A sample reference signal is shown in Fig. 3 and the normalized reference signal for each force clamp level is shown in Fig. S1 of the Supplemental Material.

3.2.3 Plant

The plant comprises a series combination of the length actuator, muscle tissue and force transducer. The plant is a single-input, single-output system, where the input is length in volts and the output is force in volts. The muscle tissue is assumed to be the main contributor to the dominant dynamics of the plant, particularly during the phase of the force clamp in which the tissue’s shortening velocity is measured since, in this phase, the effects of resonance of the force transducer and length controller in response to the initial rapid length change and the ensuing interaction with the feedback controller have subsided. This is why, in Section 3.3, modelling efforts of the plant are focused on the muscle tissue itself.

3.2.4 Feedback Controller

Feedback control helps maintain good tracking performance in the presence of disturbances and model uncertainty [20, 21]. The feedback controller is a PI controller, for which the control effort $u_{fb}(t)$ is the sum of two terms

u_{fb}(t)=k_{p}e(t)+k_{i}\int_{0}^{t}e(\tau)d\tau,

(3)

where $k_{p}$ is the proportional gain, $k_{i}$ is the integral gain and $e(t)$ is given by (1) [21]. PI controllers are ideal for tracking step reference signals because they lead to zero steady-state error thanks to the integral action [21]. The step nature of the force clamp reference signal is therefore the main reason for choosing the PI controller as the feedback controller. Prior experiments focused on finding acceptable gains showed good settling time and overshoot results with $k_{p}$ $=$ $0$ and $k_{i}$ $=$ $0.04$ s^-1.

3.2.5 Feedforward Controller

In theory, if a model perfectly describes the input-output relationship of a plant and the initial conditions are perfectly known, an inversion-based feedforward controller, which at its simplest is the inverse model of the plant, perfectly tracks the reference signal without the need for a feedback controller. However, no such model exists in practice and a feedback controller is also necessary. Knowing that the PI controller is linear and the plant is nonlinear suggests that a nonlinear feedforward controller could realize the best performance, but this will depend on the degree of nonlinearity of the plant. Using nonlinear feedforward with linear feedback is a standard approach in the field of robotics [21]. The feedforward controller is generated using a model of the plant dynamics that is estimated in real time via experimental input-output data and system ID techniques (Section 3.3.5).

3.3 System ID of Input-Output Data

3.3.1 Models

As the rapid dynamics of muscle tissue pose the greatest control challenge, modelling efforts for the feedforward controller focused predominantly on reliably capturing these dynamics. It is assumed that the contractile dynamics of the tissue change slowly relative to the time scale of a typical force clamp, allowing the tissue to be modelled as a pseudo-static, viscoelastic material, where actively changing tissue properties can be assumed to be constant over the course of a typical force clamp. The aim is to find a model that optimizes the ability to capture the muscle tissue’s properties, while keeping the number of free parameters low. This maximizes feedforward effort thus minimizing feedback effort, while also ensuring robustness and repeatability. The choice of model for the feedforward design can be reduced to the choice between linear and nonlinear models. A generalized linear model and a minimalistic nonlinear model were used to assess whether the nonlinearity is essential to achieving the control objectives.

The first model, which will be referred to as the linear model, is a generalized linear model in the form of a biproper transfer function of order $n$ given by

H(z)=\dfrac{\alpha_{n}+\alpha_{n-1}z^{-1}+\ldots+\alpha_{1}z^{-n+1}+\alpha_{0}z^{-n}}{1+\beta_{n-1}z^{-1}+\ldots+\beta_{1}z^{-n+1}+\beta_{0}z^{-n}},

(4)

where $\alpha_{i}$ and $\beta_{i}$ are the parameters to be estimated and the initial conditions are assumed to be quiescent (refer to Section 2.2.1 of the Supplemental Material for the derivation) [22]. This model structure was informed by two linear, viscoelastic models that were previously explored: the generalized Maxwell model [23] and the Golla-Hughes-McTavish model [24], both of which yield biproper transfer functions. Systems that can be represented by biproper transfer functions are known to exhibit feedthrough, i.e., a finite, non-zero gain can be observed at infinitely high input frequencies. In the case of muscle tissue, feedthrough is observed as a step change in length results in an immediate, finite change in force (see Fig. 4).

The second model, which will be referred to as the nonlinear model, is a first-order Maxwell model [25] with a zeroth-order, nonlinear spring (Fig. 5). The zeroth-order spring is chosen to be nonlinear because it is responsible for the steady-state force value following a step change in length. More specifically, it takes the form of a power law because a power law relationship was observed when plotting elastic recoil length versus force clamp values for a given muscle tissue using data from [14]. The ordinary differential equation of this model is given by

\begin{split}\dot{F}(t)=(k_{2}+k_{1}n(x(t)+L_{0})^{n-1})\dot{x}(t)-\dfrac{k_{2}}{c}F(t)\\ +\dfrac{k_{1}k_{2}}{c}(x(t)+L_{o})^{n},\end{split}

(5)

where spring constants $k_{1}$ and $k_{2}$ , viscous damping coefficient $c$ and exponent $n$ are to be estimated (refer to Section 2.2.2 of the Supplemental Material for the derivation). The muscle tissue’s projected unstrained length $L_{o}$ is given by

L_{o}=\left(\dfrac{F(0)}{k_{1}}\right)^{\!1/n}.

(6)

3.3.2 Designing Input-Output Data for Model Estimation and Validation

Separate datasets are used for model estimation and validation as dynamic tissue remodelling can cause past length changes to affect future force responses. The input-output data need to capture the dynamics relevant to the control of force clamps, which combine a rapid step-like phase with a maintained, slower phase. To capture these dynamics, the input signal is designed using ramped step length changes of a range of amplitudes. The upper bound of the input is set to 0 V, which corresponds to the length at which the tissue is isometrically contracted, i.e., the reference length, such that the muscle tissue is never stretched beyond this length to avoid tissue damage. The lower bound of the input is set to the upper bound minus 2 $\%$ of the muscle tissue’s reference length to avoid slack muscle tissue which is problematic for system ID. Sample input-output data to estimate the nonlinear model are shown in Fig. 4. Refer to Figs. S2 and S3 of the Supplemental Material for sample input-output data to estimate the linear model and sample input-output data to validate both models, respectively.

The input-output data are preprocessed prior to model estimation and validation in order to yield the best results for a given model. There is a delay of approximately 1 ms between the input and output data caused by response delays of the length actuator and force transducer, which cannot be captured by the models. The delay is removed by shifting the output data relative to the input data for both model estimation and validation. For the linear model, since (4) is formed on the assumption that the initial conditions are quiescent, the offset is removed from the output data for model estimation and validation so that the starting force is 0 V. Also, in order for the computation of the least squares estimate (Section 3.3.3) to be numerically well-conditioned, the input-output data for model estimation are normalized to the range $\left[-1,1\right]$ . Prior to model validation, the scaling is reversed to obtain the original transfer function $H(z)$ . For the nonlinear model, the output data for model estimation and validation are filtered using a moving average filter with a window size of 20 samples in order to smooth out any resonance and measurement noise from the force transducer.

3.3.3 Model Estimation

The parameters of the linear model are estimated for order $n\leq 3$ using MATLAB’s tfest [26], which involves least squares in its default settings. Although higher order models have the potential to improve the quality of fit, they give rise to more variability in the parameter estimates as a result of overfitting. The resulting discrete-time transfer function is then converted to continuous time by tfest.

The parameters of the nonlinear model are estimated via constrained optimization using MATLAB’s fmincon [27] in its default settings. Equation (5) is solved at every iteration of fmincon to minimize the cost function given by the normalized root mean square error (NRMSE),

\textrm{NRMSE}=\sqrt{\dfrac{\dfrac{1}{N}\sum_{i=0}^{N}(y_{i}-\hat{y}_{i})^{2}}{\dfrac{1}{N}\sum_{i=0}^{N}y_{i}^{2}}},

(7)

where $\boldsymbol{y}$ is the observed output, $\hat{\boldsymbol{y}}$ is the simulated output and the error or the residual $\boldsymbol{r}$ is defined as $\boldsymbol{r}=\boldsymbol{y}-\hat{\boldsymbol{y}}$ . The NRMSE is equal to 0 when there is complete overlap between the simulated and observed outputs. In this case, the observed output is the preprocessed observed output. Equation (5) is solved using the classical fourth-order Runge-Kutta method [28] since a fixed time step is needed for fmincon to converge. The initial guess and constraints for each of the free parameters of the nonlinear model are provided in Section 2.2.2 of the Supplemental Material.

3.3.4 Model Validation

Regardless of the type of models used for this application, NRMSE is used to assess 1) the quality of fit between the simulated output and the preprocessed observed output for model validation and 2) the quality of fit between the total control effort and the feedforward signal for each force clamp level. For the linear model, NRMSE is also used to select the model order that results in the best fit, i.e., the model order for which the NRMSE is closest to 0 using the simulated output and the preprocessed observed output for model validation.

3.3.5 Feedforward Controller

Once the continuous-time transfer function of the linear model is obtained, it is inverted by switching the numerator and denominator in order to generate the transfer function of the feedforward controller. However, if prior to the inversion the numerator contains nonminimum phase zeros, which are zeros in the open right half plane, they are first mirrored about the imaginary axis into the open left half plane [20]. This is done to ensure that the transfer function is bounded-input, bounded-output stable once inverted, i.e., all of its poles are in the open left half plane. Mirroring non-minimum phase zeros does not change the gain of the transfer function, but it changes its phase. The feedforward signal is then simulated using MATLAB’s lsim [29] with the reference signal as the input to the feedforward controller. A median filter is applied to remove an initial spike in the response due to the drop in force since the linear model does not accurately estimate the feedthrough.

Once the parameters of the nonlinear model are estimated, they are substituted into the ordinary differential equation of the feedforward controller given by

\begin{split}\dot{x}(t)=\dfrac{1}{(k_{2}+k_{1}n(x(t)+L_{o})^{n-1})}(\dot{F}(t)+\dfrac{k_{2}}{c}F(t)\\ -\dfrac{k_{1}k_{2}}{c}(x(t)+L_{o})^{n}),\end{split}

(8)

which is the inverse input-output relationship of (5) (refer to Section 2.2.2 of the Supplemental Material for the derivation). The feedforward signal is then simulated using the reference signal as the input to the feedforward controller.

3.4 Data Analysis

A sampling frequency of 10 kHz is used, which is more than the minimum required of ten times the assumed bandwidth [30] of muscle tissue [19]. Prior to calculating the settling time, the force response is normalized by the reference contractile force. It is also filtered using a moving average filter with a window size of 20 samples; otherwise, the measurement noise affects the settling time calculation. As shown in Fig. 3, the settling time is calculated as the time it takes for the force response to settle within an error band from the start of the force clamp at $t_{0}~{}=$ 0.08 s. An error band of $\pm$ 0.25 $\%$ of the starting force is chosen because it corresponds to less than a $\pm$ 20 $\%$ deviation of the measured shortening velocity for a 5 $\%$ force clamp, which is acceptable for the research needs in this study. A constant absolute magnitude error band is used as control effort is proportional to the absolute error and for higher force clamps, measurement noise would exceed a constant relative magnitude error band. Using Fig. 3 as an example, the overshoot $OS$ is calculated as

OS=\dfrac{(1-\frac{F}{F_{\text{ref}}})-(1-\frac{F_{\text{level}}}{100})}{1-\frac{F_{\text{level}}}{100}}\times 100\%,

(9)

where $F_{\text{level}}$ is the force clamp level. Note that an overshoot of 5 $\%$ for a 5 $\%$ force clamp corresponds to slack muscle tissue.

4 Results

4.1 System ID Results

Comparison of the linear versus the nonlinear models, using the same input-output data, clearly shows better predictive power of the nonlinear model in open loop for each data point (Fig. 6 (a)); sample output data for model validation is provided in Fig. S4 of the Supplemental Material. This is further confirmed when comparing the quality of the feedforward signal generated from both models relative to the total control effort from the 2DOF control configuration, which uses the feedforward signal generated by the nonlinear model. Fig. 6 (b) shows the NRMSE of this comparison for the 5 $\%$ force clamps, while Fig. S5 of the Supplemental Material shows a sample of the comparison for 5 $\%$ and 80 $\%$ force clamps. The larger the amplitude of the force clamp level, the greater the display of the muscle tissue’s nonlinear, viscoelastic properties, which the linear model falls short in capturing as compared to the nonlinear model. These data clearly show that the nonlinear model performs better than the linear model in reducing the control effort of the PI controller. To confirm the validity of the force clamps that correspond to the total control efforts used to compute NRMSE, their settling time and overshoot values are provided in Table S1 of the Supplemental Material. From this point forward, the remaining results will only be presented with respect to the nonlinear model.

The repeatability of the estimation process for the nonlinear model was assessed by estimating this model three consecutive times within the same contraction for ASM (Fig. 7 (a)) and across contractions with the same reference contractile force for FDB muscle (Fig. 7 (b)). If the estimation process yields repeatable results for a tissue with a given reference contractile force, then the feedforward signal for this tissue for a given force clamp should be fairly consistent. However, it can be observed from the average NRMSE values presented in these figures that there is some variability. The corresponding feedforward signals of the ASM tissue with the greatest variance based on the 95 $\%$ confidence intervals in Fig. 7 (a) are displayed in Fig. 7 (c). The variability can be observed in the parameter estimates within and across contractions for all muscle tissues and it manifests itself in the feedforward signals that are generated. The coefficient of variation for the parameters of each model estimate in Fig. 7 (a) is shown in Fig. 7 (d). Note that the coefficient of variation can get distorted for small sample means as is the case for $c$ , $k_{2}$ and sometimes $k_{1}$ . For either ASM or FDB muscle, the greatest variability is observed for $k_{1}$ and $c$ .

The quality of the feedforward signal generated by the nonlinear model for repeated force clamps at the same contractile force and at reduced contractile force in FDB muscle is shown in Fig. 8. This was done in FDB muscle only because it is easier to reliably reduce the contractile force with EFS contractions. The NRMSE values in Fig. 8 are computed using the reference signal and force output for each 5 $\%$ force clamp. Lower NRMSE values are achieved when the feedforward signal is generated using the model that is estimated at the same contractile force as the force clamp, which suggests that the system ID process may need to be repeated for large differences in contractile force. Table S2 of the Supplemental Material, which includes the corresponding values of settling time and overshoot, supports this observation and also suggests that a change in control gains may be required in addition to model re-estimation.

4.2 Force Control Results

The average settling time and average overshoot are shown for repeated 10 $\%$ , 7 $\%$ and 5 $\%$ force clamps for ASM and FDB muscle in Fig. 9. A trade-off between settling time and overshoot is observed, especially in Fig. 9 (c), where, for example, the faster the settling time, the greater the overshoot. A user could adjust the control gains depending on their overshoot versus settling time requirements. In the majority of cases, the control objectives regarding these two parameters are met. The excessive settling times corresponded to a high NRMSE in validation, i.e., higher model uncertainty (Fig. 6 (a)). The inter-tissue variability that is observed for either muscle type suggests that some aspect of the muscle dynamics is not fully captured, yet variable between tissues. This may be countered by manually adjusting $k_{i}$ to meet given research requirements for settling time and overshoot, which was done here for two FDB muscle tissues in order to maintain closed-loop stability ( $k_{i}$ was lowered from 0.04 s^-1 to 0.02 s^-1). Refer to Figs. S6 and S7 of the Supplemental Material for sample force responses for each force clamp level and muscle tissue type.

To assess the contribution of the feedforward controller, the 2DOF control configuration is compared to the one-degree-of-freedom control configuration with only a feedback controller for the 10 $\%$ , 7 $\%$ and 5 $\%$ force clamp levels in Fig. 10. Regardless of the type of muscle tissue, the 2DOF control configuration results in faster settling times. Even though the overshoot is higher for the 2DOF control configuration, the overshoot remains below the set limit of 3 $\%$ and the priority is in reducing the settling time.

5 Discussion

Modelling efforts were geared towards generating a nonlinear inversion-based feedforward controller to better track the reference signal for each force clamp level, reducing or eliminating the need for supervision from the user. The nonlinear model captures the nonlinear, viscoelastic properties of the muscle tissues better than the linear model. However, the results indicate that improvements can still be made in terms of capturing these properties and ensuring the repeatability of the model estimation. It is currently possible to use the same model estimate for force clamps at the same contractile force, but with reduced quality for force clamps at different contractile forces, which supports the assumption that a given model estimate is only valid for a given pseudo-static state. As system ID may need to be repeated whenever the contractile force is substantially changed, supervision may be required from the user to judge whether or not a new model estimate is needed.

As shown in Fig. 7 (d), there can exist a large degree of variability when estimating model parameters. Imposing tighter constraints on either the nonlinear parallel stiffness parameter $k_{1}$ or the viscosity parameter $c$ as well as exploring different cost functions for model estimation could reduce the variability that is observed in the parameter estimates within and/or across contractions for all muscle tissues. One could also try increasing the maximum amplitude of the estimation and validation inputs for system ID to ensure that the system ID process captures more of the nonlinear behaviour of the force-length relationship of the muscle tissue. Because of the shortening velocity’s sensitivity to excessive overshoot, the feedforward signal should not overestimate the total control effort for any force clamp level as the integral-only design of the feedback controller does not react fast enough to reduce overshoot caused by an overestimating feedforward signal.

One set of improvements to the control was tested to see if fixing the value of $c$ could improve repeatability. This was combined with the introduction of a regularization term of the form $\alpha\lVert\mbf{p}\rVert_{2}^{2}=\alpha\mbf{p}^{\top}\mbf{p}$ , where $\alpha$ is a scaling factor and $\mbf{p}=[n$ $k_{1}$ $k_{2}]^{\top}$ , using the same cost function in order to reduce the impact of the trade-off between high bias and high variance when estimating the model parameters. In comparing Figs. 7 and 11, one can see that the variability in the parameter estimates is greatly reduced with these improvements in place.

The current force control scheme is robust for different tissue types, as shown by the similarities in the results for FDB muscle and ASM. Moreover, improving the repeatability of the system ID process, specifically that of model estimation, should in turn improve the repeatability of the force control scheme. In regards to the gains of the PI controller, there are methods to automatically tune them with minimal manual tuning required, e.g., Ziegler-Nichols tuning method in [21]. However, for this application, given that the control objectives specific to settling time and overshoot are met for most cases, some minor improvements in the system ID process as suggested above may be sufficient to allow just one fixed set of control gains. Nonetheless, if the users have the option to fine-tune the control gains, they would only need to fine-tune $k_{i}$ given $k_{p}$ $=0$ .

6 Conclusion

This paper introduces the use of 2DOF control with nonlinear feedforward through system ID in the force control of biological tissues and the approach can be generalized to any force control of nonlinear, viscoelastic materials. The force control scheme comprises a feedback PI controller and an inversion-based feedforward controller leveraging a nonlinear physiological model in a system ID context that is superior to classic linear system ID. Results from experiments with equine ASM and murine FDB muscle are presented to assess the performance, robustness and repeatability of the force control scheme as well as the number of inputs and level of supervision required from the user. While improvements can still be made as part of future work for this application, the force control scheme was able to fulfill the stated control objectives in most cases, including the requirements for settling time and overshoot. In addition, it paves the way for using force control in long-term incubation studies because of the greater level of autonomy that it offers.

Acknowledgment

This work was supported by the Natural Sciences and Engineering Research Council of Canada and Aurora Scientific Inc. through a Collaborative Research and Development grant for the Lauzon Lab, and by the Fonds de recherche - Nature et technologies Master’s Research Scholarship awarded to Amanda Bianco.

References

[1] S. Munevar et al., “Traction force microscopy of migrating normal and h-ras transformed 3t3 fibroblasts,” Biophysical journal, vol. 80, no. 4, pp. 1744–57, 2001.
[2] I. M. Tolić-Nørrelykke et al., “Spatial and temporal traction response in human airway smooth muscle cells,” American journal of physiology. Cell physiology, vol. 283, no. 4, pp. 1254–66, 2002.
[3] H. Matsumoto et al., “Comparison of gel contraction mediated by airway smooth muscle cells from patients with and without asthma,” Thorax, vol. 62, no. 10, p. 848, 2007.
[4] J. Chamley-Campbell et al., “The smooth muscle cell in culture,” Physiological Reviews, vol. 59, no. 1, pp. 1–61, 1979.
[5] H. A. Arafat et al., “Heterogeneity of bladder myocytes in vitro: modulation of myosin isoform expression,” Tissue and Cell, vol. 33, no. 3, pp. 219–232, 2001.
[6] A. Huber and S. F. Badylak, “Phenotypic changes in cultured smooth muscle cells: limitation or opportunity for tissue engineering of hollow organs?,” Journal of Tissue Engineering and Regenerative Medicine, vol. 6, no. 7, pp. 505–511, 2012.
[7] C. Y. Seow, Introduction to Smooth Muscle Mechanics. FriesenPress, 2018.
[8] C. Y. Seow, “Hill’s equation of muscle performance and its hidden insight on molecular mechanisms,” The Journal of general physiology, vol. 142, no. 6, pp. 561–73, 2013.
[9] T. A. McMahon, Muscles, reflexes, and locomotion. Princeton paperbacks, Princeton, N.J.: Princeton University Press, 1984.
[10] O. S. Matusovsky et al., “Peripheral airway smooth muscle, but not the trachealis, is hypercontractile in an equine model of asthma,” American journal of respiratory cell and molecular biology, vol. 54, no. 5, pp. 718–27, 2016.
[11] O. S. Matusovsky et al., “Cd4+ t cells enhance the unloaded shortening velocity of airway smooth muscle by altering the contractile protein expression,” The Journal of physiology, vol. 592, no. 14, pp. 2999–3012, 2014.
[12] S. R. Bullimore et al., “Could an increase in airway smooth muscle shortening velocity cause airway hyperresponsiveness?,” American journal of physiology. Lung cellular and molecular physiology, vol. 300, no. 1, pp. 121–31, 2011.
[13] L. Luo et al., “The huxley crossbridge model as the basic mechanism for airway smooth muscle contraction,” American journal of physiology. Lung cellular and molecular physiology, vol. 317, no. 2, pp. L235–L246, 2019.
[14] G. Ijpma et al., “Intrapulmonary airway smooth muscle is hyperreactive with a distinct proteome in asthma,” The European respiratory journal, 2020.
[15] P. J. Sabourin et al., “An active feedback system for isotonic studies of smooth muscle,” IEEE Transactions on Biomedical Engineering, vol. 38, no. 6, pp. 614–616, 1991.
[16] G. Ijpma et al., “Maintenance of contractile function of isolated airway smooth muscle after cryopreservation,” American journal of physiology. Lung cellular and molecular physiology, vol. 315, no. 5, pp. L724–L733, 2018.
[17] A. R. Demonbreun and E. M. McNally, “DNA electroporation, isolation and imaging of myofibers,” Journal of visualized experiments : JoVE, no. 106, p. e53551, 2015.
[18] P. Keire et al., Isolation and Culture of Skeletal Muscle Myofibers as a Means to Analyze Satellite Cells, pp. 431–468. Totowa, NJ: Humana Press, 2013.
[19] R. F. Kirsch et al., “Muscle stiffness during transient and continuous movements of cat muscle: perturbation characteristics and physiological relevance,” IEEE Transactions on Biomedical Engineering, vol. 41, no. 8, 1994.
[20] L. Qiu and K. Zhou, Introduction to Feedback Control. Prentice Hall, 2010.
[21] K. J. Åström and T. Hägglund, Advanced PID control. Research Triangle Park, NC: ISA-The Instrumentation, Systems, and Automation Society, 2006.
[22] J. P. Hespanha, “Advanced undergraduate topics in control systems design.” https://web.ece.ucsb.edu/%7Ehespanha/published/allugtopics-20220327.pdf, 2019.
[23] D. Gutierrez-Lemini, Engineering Viscoelasticity. Boston, MA: Springer US, 2014.
[24] D. J. McTavish and P. C. Hughes, “Modeling of linear viscoelastic space structures,” ASME Transactions Journal of Vibration Acoustics, vol. 115, pp. 103–110, 1993.
[25] Y. C. Fung, Biomechanics : mechanical properties of living tissues. New York: Springer-Verlag, 2nd ed. ed., 1993.
[26] MathWorks, “tfest.” https://www.mathworks.com/help/ident/ref/tfest.html. (accessed 2019).
[27] MathWorks, “fmincon.” https://www.mathworks.com/help/optim/ug/fmincon.html. (accessed 2020).
[28] M. T. Heath, Scientific computing: an introductory survey. Boston: McGraw-Hill, 2nd ed. ed., 2002.
[29] MathWorks, “lsim.” https://se.mathworks.com/help/control/ref/lti.lsim.html. (accessed 2019).
[30] L. Ljung, System identification: theory for the user. Prentice Hall information and system sciences series, Upper Saddle River, NJ: Prentice Hall PTR, 2nd ed. ed., 1999.