Abstract
Small diameter tissue-engineered arteries improve their mechanical and functional properties when they are mechanically stimulated. Applying a suitable stress and/or strain with or without a cycle to the scaffolds and cells during the culturing process resides in our ability to generate a suitable mechanical model. Collagen gel is one of the most used scaffolds in vascular tissue engineering, mainly because it is the principal constituent of the extracellular matrix for vascular cells in human. The mechanical modeling of such a material is not a trivial task, mainly for its viscoelastic nature. Computational and experimental methods for developing a suitable model for collagen gels are of primary importance for the field. In this research, we focused on mechanical properties of collagen gels under unconfined compression. First, mechanical viscoelastic models are discussed and framed in the control system theory. Second, models are fitted using system identification. Several models are evaluated and two nonlinear models are proposed: Mooney-Rivlin inspired and Hammerstein models. The results suggest that Mooney-Rivlin and Hammerstein models succeed in describing the mechanical behavior of collagen gels for cyclic tests on scaffolds (with best fitting parameters 58.3% and 75.8%, resp.). When Akaike criterion is used, the best is the Mooney-Rivlin inspired model.
1. Introduction
In vascular tissue engineering (VTE), bioreactors are used to subdue cells in culture on scaffolds with quasi-physiological conditions during the maturation process that is expected to induce extracellular matrix reorganization and to provide mechanical properties to the regenerated tissue thus leading to the development of functional tissue [1, 2]. Strains can easily be generated and measured in bioreactors but stresses can only be estimated, and this is not a trivial task. Commercial bioreactors can generate and noninvasively measure pressure and diameter of constructs but measurements of stresses can only be estimated by models. One of the most promising approaches for the estimation of the stresses constitutes in generating a mechanical model [3]. Mechanical models are needed not only to identify the degree to which the regenerated tissue will match the physiological tissue [4], but also to provide an estimation of the generated internal stress and its role in growth and development of the constructs [1]. Scaffolds made of synthetic or natural biodegradable materials are generally used, and there is a considerable interest in tailoring mechanical properties of scaffolds to facilitate the cell growth [5]. Among the natural scaffolds, collagen gels have received special attention mainly because they present unmatched biological performances, such as its nontoxicity, low immunogenicity, and antigenicity [6]. In addition to its biological relevance, the collagen is the most abundant protein in native tissues and can be isolated easily in large quantities. Specifically, the main functional feature of collagen is load bearing of tensile force in most tissues where mechanical function is essential such as blood vessels, skin, cartilage, and bones [7, 8]. The development of models of the mechanical behaviour of collagen will also facilitate the optimization of the mechanical properties of the construct by modifying the conditions generated in the bioreactor [4].
Stress and strain present a nonlinear behaviour in most of the biological and bioartificial tissues [9–11]. In particular, viscoelasticity has an important role in these materials, mainly because the entire history of the straining determines the material properties. In addition, there is a dependence of the loading rate in the mechanical response. In the search of mathematical relations between force and displacement Fung proposed the quasilinear viscoelastic model [11]. Nekouzadeh et al. [12] and Nekouzadeh and Genin [4] described the two main limitations of this approach: the accuracy and the computational cost in calibration experiments. Therefore the selection of the nonlinear viscoelastic model and the fitting procedure for a specific tissue can also be a computational challenge.
Although the mechanical model of other scaffolds used in VTE (e.g., standard polyglycolic acid gels) is rather well known [2, 13], for collagen gels, the model describing the mechanical behaviour has much less been studied. Shear and uniaxial extension experiments have been studied by Meghezi et al. [10], using relaxation tests, thus proposing a viscoelastic model, represented by a spring (elastic modulus) and Maxwell elements, associated in parallel. A generalization of Fung’s model based on incremental uniaxial tests (step and ramp) where both the spring moduli and time constants vary with strain was also proposed by Pryse et al. [14]. Chandran and Barocas have proposed a review of the mechanical properties of collagen gels [9]. The main conclusion was that the collagen network is actively involved in shear and tensile tests. On the contrary, for compression the mechanisms are not well established yet. Achilli and Mantovani [5] measured the mechanical properties (Young modulus) under confined-compression in a study aimed to gel characterization after preconditioning and assuming a linear model. Chandran and Barocas [9] observed that under step-confined compression (10%), collagen gels exhibited a collapse while, under ramp-compression (0.1%/second), no collapse was observed. They concluded that the gel behaves primarily as a viscoelastic solid, with important damping but negligible relaxation on the time scale of their experiments (1800–2400 seconds).
This work focuses on the mechanical constitutive equations of collagen gels in cyclic compression (loading and unloading). We assumed that the mechanical properties of collagen gels should be evaluated under cyclic testing in order to better mimic the cyclic solicitations imposed by the bioreactor to the vascular construct [15]. Mechanical models are written in the form of a dynamic system to frame the work within the systems control theory [16]. In control engineering, the field of system identification (SI) is devoted to building mathematical models of dynamical systems from measured data [17]. Using SI, we implemented a simple and fast way of achieving constitutive equations describing the behaviour of collagen gels under unconfined compression tests. Results were then interpreted by linear viscoelastic models [11], the Moolin-Rivlin approach [18], and a new polynomial model is finally proposed to describe the mechanical behaviour of collagen gels.
2. Material and Methods
2.1. Sample Preparation
The protocol used to prepare the gels was previously described in [19]. Briefly, type I collagen is extracted from rat-tail tendons and solubilized in acetic acid solution (0.02 N) at 4 g/L. A collagen solution (2 g/L) is mixed with Dulbecco’s modified Eagle medium (DMEM, Gibco, Invitrogen Corporation, Burlington, ON, Canada, 1.1X), NaOH (15 mM), and HEPES (20 mM) in deionized water. This mixture is poured into moulds and left to jellify overnight at 4°C.
2.2. Mechanical Tests
Disk-shaped samples were tested in unconfined compression mode. The 70 mm diameter Teflon compression plate is attached to a 10 N load cell mounted on an Instron Microtester (Instron 5848 Microtester, Instron Corporation, Norwood, MA, USA) as shown in Figure 1. The unconfined compression was performed with a uniform ramp speed ranging from 0.1 to 1 mm/min. The average strain was 20%. Eight cycles (15 to 25% strain) were performed followed by a relaxation test. For this final test, the strain rate is increased to 10 mm/min with a final strain value of 50%. Relaxation tests with final strains of 25% and 50% were also performed with strain rate of 10 mm/min. All tests were performed at 37°C.
2.3. Mechanical Models
In this section we will focus on the two proposed mechanical models. Basically, they consist in two Maxwell bodies [11] with nonlinear springs. Two nonlinear springs relationships will be tested: Mooney-Rivlin and polynomial (models 1 and 2, resp.). The general approach is quasi-static, and acceleration effects are neglected. Consequently, the inertia parameter is neglected in the evaluated models.
2.3.1. Model 1
The first model is inspired in the Mooney-Rivlin model, which is applied to uniaxial compression of incompressible and isotropic gels. It leads to a stress-stretch relation written aswhere is the Eulerian stress and is the stretch ratio ( is the extension and is the extension at instant 0). and are constants. Equation (1) yields a linear relationship between and . Note also that this becomes the Neo Hookean model when [18].
The two Maxwell bodies proposed (with the modified Mooney-Rivlin model) have the next equations:where and are the extension of the dashpot and the spring , respectively. The extension is defined as is the spring constant and the viscosity. Finally, the load is written as follows:where and is the cross-section area of the sample. It should be remarked that the extension varies with time but does not depend on the derivatives of . By definition, it is a static function of thus belonging to a memoryless system. The stretch is also a static function of .
2.3.2. Model 2
A linear representation of two Maxwell bodies iswhere is the strain, and is the relaxation time due to Maxwell body “”. Instead of using , we propose here to use a nonlinear function of the strain, in particular, an order three polynomial function of the strain:
2.4. Nonlinear System Modeling
In this section we will consider the above described viscoelastic models as nonlinear dynamical systems in state space representation. From now on we will denominate , , and as input, output, and state space variable, respectively. All the mentioned variables vary with time.
2.4.1. Model 1
Replacing in (2)–(4) the state space representation variables ( and ) and taking and , a nonlinear system is now defined aswhere expressed in function of the input of the system is
2.4.2. Model 2
On the other hand, for implementing the polynomial model we define and , and (5) is modified towhere , , , and . The state space representation variables are now defined as and , and a state space system is obtained (see [16]):where and are the elements of the row vector . This is a linear system but if instead of using we use , a block-oriented Hammerstein nonlinear system is obtained (see Figure 2).
Equations (7) and (10) (using as input) constitute the nonlinear systems to be tested in the following sections.
2.5. System Identification
The fitting procedure of (7) and (10) to real data is a computational challenge. SI is a well established area of control system theory devoted to find mathematical models from measurements [17, 20, 21]. The basic idea of these techniques is as follows. Let us assume that an input signal is applied to a dynamical system for a determined time interval . The corresponding output signal is obtained, . Suppose that both signals are sampled at instants, then data can be denoted by . A general form (or structure) for the mathematical model is assumed and then the parameters are determined from experimental data. Often, a variety of model structures are evaluated, and the most successful ones are retained. The data are used to build the model by minimizing the difference between the output of the estimated model and the measured output . This minimization can be based on statistical techniques which ensure the robustness of the method. An important step in a SI process is the validation of the model, which must be performed with an independent data set. In this work, we will use one sample gel data set for identification of the model and the remaining samples for validation purposes. One of the validation parameter iswhere and are the output given by the estimated model and the measured data, respectively. It is important to remark that the validation is not a new fitting: the measured input is applied to the identified model and then the obtained output is compared to the measured one. The final prediction error (FPE) is also defined in Appendices as another parameter for evaluation of the model. This parameter is computed with the identification data.
Several tools were developed to implement the SI algorithms. In this paper two Matlab toolboxes are used: System Identification and CON-tinuous-Time System IDentification (CONTSID toolbox [22]).
Particularly in this work, linear and nonlinear SI are tested. For the former, srivc algorithm is applied [20] using the structure given by (10) and the strain as input (). For the nonlinear SI we implement grey-box SI using (7) (model 1) and (10) (model 2, using (6) as input ). Levenberg-Marquardt algorithm was used to minimize the prediction error of the models.
3. Results
3.1. Relaxation
Figure 3 shows the relaxation tests performed at two strain steps, 25% (one sample) and 50% (three samples). Continuous time linear system identification, particularly, srivc algorithm was used to obtain transfer functions of the form of (A.1). The input of the system is the step of the strain () and the output is the load (). The best model outputs are also plotted in Figure 3. Curve almost overlaps for the 25% test and the errors (equal to which are represented in the inset figure) were always lower than 0.01 for the 50% tests. The dynamics of such identified model and the physical Maxwell bodies time constants () are related. In fact, the poles of the transfer functions are the inverse of the time constants of the Maxwell bodies. The values obtained for 25% of strain are and seconds and for 50%, and seconds. The results presented here give rise to the proposal of the two Maxwell bodies model for this study.
(a)
(b)
3.2. Cycling
Table 1 shows the experimental conditions of the cyclic tests for the four samples. As indicated by (1) in Section 2.3.1, a plot of versus helps in determining whether Mooney-Rivlin model fits some region of the stress-strain behavior. Such a plot in Figure 4 suggests that the first compression (plotted in dashed-dotted black line) can be represented by a Neo Hookean model (which is plotted in dashed blue line). The compression phases of the subsequent cycles (plotted in red line) are well described by Mooney-Rivlin model (plotted in dotted blue line). The hysteresis loop appears because of the cyclic test and is an evidence of a viscous response of collagen gels. These findings give rise to the proposal of the Mooney-Rivlin inspired model described by (7).
Identification and validation data are presented in Figure 5. In order to excite all the dynamical modes, the fastest strain data are used for identification which correspond to sample 1.
(a)
(b)
(c)
(d)
The best linear model identified is presented in Figure 6. The model structure is given by (A.1). (strain) is used as input and as output. The continuous time srivc algorithm is used. The parameters obtained are s and s.
(a)
(b)
(c)
Figure 7 shows the results for the Mooney-Rivlin model. Algorithms time constants were the ones identified in the linear cases ( s and s). For model 1, the shape of the wave is followed almost superimposed and the estimated parameters are , , , and .
(a)
(b)
(c)
If the constant is now fixed to 0, the model 1 is now similar to a Neo Hookean model. Figure 8 shows the results when the identification of the model is performed with such a model structure. The first compression of sample 1 is well reproduced by the model 1 (with ). This is not the case of the validation data Figures 8(b) and 8(c). The estimated parameters are , , and .
(a)
(b)
(c)
The results for the Hammerstein model are presented in Figure 9. The parameters obtained were and 0 for , and 2, respectively, and also , , s, and s.
(a)
(b)
(c)
The FIT parameters obtained for the best models are presented in Table 2 for the validation data. The Hammerstein model yields the highest FIT parameter: 75.8%. It occurred for sample 3. For sample 4, Hammerstein model also yields a high FIT parameter. Finally, the FPE parameters (computed with identification data at cycling part only) for the linear model, model 1 (Mooney-Rivlin inspired model), and model 2 (Hammerstein with polynomial function of the strain) are , , and , respectively.
4. Discussions
This study explored several aspects of the nonlinear relation between uniaxial extension and load for collagen based gels. Compression relaxation and cyclic tests were realised in order to study the mechanical behavior. System identification techniques for processing data were proposed.
Relaxation compression tests indicate that two Maxwell bodies can represent the collagen gel presented here. Regarding the cyclic response of the material, linear Maxwell behavior can not correctly explain the measurements. To interpret it, the analysis of the cyclic tests is divided in two parts: (i) first excursion of the cyclic test (first compression) and (ii) the rest of the cycling.
Figure 8 shows the results of the first excursion (i). According to the identification data, the Neo Hookean inspired model (model 1 using ) seems to represent very well the measurements (Figure 8(a)), but this must be carefully interpreted. It should be noted that the fitting of validation data for the first excursion is not good (Figures 8(b) and 8(c)). Strictly, it can not be concluded that this is the mechanical model of the gels.
For the rest of the cycling data (ii), the linear model (Figure 6) do not represent qualitatively the measurements. However the FIT parameters are surprisingly high. The Mooney-Rivlin inspired model (model 1) has a relatively good FIT parameter, 58.3% for sample 3 (Table 2). If we observe only the second part of the cycling (from 200 to 600 seconds), the fitting is still better (curves almost overlap for samples 2 and 3, see Figure 7(b)). This is not the case of the sample 4; see Figure 7(c). With the best Hammerstein model from 200 to 600 seconds the fitting is even better (see Figure 9(b)). The Hammerstein model achieved the best fitting when it was applied to sample 4: 58.2% (Figure 9(c)). Although Mooney-Rivlin models have very good fitting parameters for this part of the compression test, Hammerstein structure, both quantitatively and qualitatively, yields a better description of the wave form. Considering the FPE parameter the best result is obtained with the Mooney-Rivlin inspired model. This indicates that according to the Akaike criterion the best result is obtained with the model 1.
A commentary should be added regarding the mechanical parameters ( and ) obtained with the presented technique. Although modified two Maxwell bodies model is the proposed model, we cannot strictly say that is the viscosity and is the modulus (note that for model 1 the extension depends on and ). However, we can compute , as samples have ; then it yields the value 3.4 kPa s. Using the same reasoning for 570 kPa s. Both values are comparable with the one obtained in [10].
Arteries are complex hierarchical systems, in which three types of cells (endothelial, smooth muscular, and fibroblasts) cohabit in perfect symphony in a sort of basal membrane composed by several types of proteins and biochemical factors all devoted to providing adequate support and specific signals for regulating the overall biological activity and, as a consequence, the mechanical properties [23]. When bioengineers target the regeneration of an artery, by definition, they focused on an engineering approach thus aiming at reproducing a pseudoartery. Cells are often the missed components, although their importance is crucial for their ability to regenerate tissue, organize proteins, and secrete their own extracellular matrix. In this context, computational model can be very effective in providing bioengineers the important factors to consider during experimental approaches.
5. Conclusions
The viscoelastic nature of collagen gels under dynamic solicitation is encountered in bioreactors or in mechanobiology experiments. This behavior is particularly hard to anticipate. The tools and models commented on here frame the discussion in terms of system identification methods. This approach introduces engineering concepts, such as the design of experiments and optimal experimental design, to increase the amount of information that can be retrieved from experiments. In future work we plan to study living cells within the collagen gel matrix and to monitor evolution of the mechanical properties using the models presented here.
Appendices
A. Transfer Function Model
Equation (9) can be transformed into a linear system representation called transfer function. It is obtained by taking the Laplace transform as follows:where now is the Laplace frequency and and are the Laplace transform of and , respectively. From (A.1) it can be seen that the dynamic of this system is dominated by the two poles ( and ) which are the inverse of and .
B. Validation Parameter
In order to evaluate the estimated models, FIT parameter was calculated. Another parameter is also computed, the Akaike final prediction error (FPE) defined aswhere is the number of estimated parameters and is the number of estimation (identification) data and is called loss function: The computation is made with the identification data.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank Richard Gagnon, Sébastien Meghezi, Lucie Lévesque, Nina Bono, and Ranna Tolouei for their kind assistance in revising the paper. This work was partially supported by Agence Universitaire de la Francofonie (AUF). The authors also acknowledge the Natural Science and Engineering Research Council of Canada (D. Mantovani), Centre d’Études et de Recherche en Matériaux Avancés (B. Drouin and D. Mantovani), and Emerging Leaders in the Americas Program for scholarship to R. M. Irastorza.