Abstract
Skeletal muscle system has nonlinear dynamics and subject-specific characteristics. Thus, it is essential to identify the unknown parameters from noisy biomedical signals to improve the modeling accuracy in neuroprosthetic control. The objective of this work is to develop an experimental identification method for subject-specific biomechanical parameters of a physiological muscle model which can be employed to predict the nonlinear force properties of stimulated muscle. Our previously proposed muscle model, which can describe multiscale physiological system based on the Hill and Huxley models, was used for the identification. The identification protocols were performed on two rabbit experiments, where the medial gastrocnemius was attached to a motorized lever system to record the force by the nerve stimulation. The muscle model was identified using nonlinear Kalman filters: sigma-point and extended Kalman filter. The identified model was evaluated by comparison with experimental measurements in the cross-validation manner. The feasibility could be demonstrated by comparison between the estimated parameter and the measured value. The estimates with SPKF showed 5.7% and 2.9% error in each experiment with 7 different initial conditions. It reveals that SPKF has great advantage especially for the identification of multiscale muscle model which accounts for the high nonlinearity and discontinuous states between muscle contraction and relaxation process.
1. Introduction
1.1. Muscle Modeling and FES
Functional electrical stimulation (FES) is an effective technique to evoke artificial contractions of paralyzed skeletal muscles. It has been employed as a general method in modern rehabilitation to partially restore motor function for patients with upper neural lesions [1, 2]. Recently, the rapid progress in microprocessor technology provided the means for computer-controlled FES systems [3–5]. A fundamental problem concerning FES is how to handle the high complexity and nonlinearity of the neuromusculoskeletal system [6, 7]. In addition, there is a large variety of patient situations depending on the type of neurological disorder. To improve the performance of motor neuroprosthetics beyond the current limited use of such system, subject-specific modeling is essential. The use of a mathematical model can improve the development of neuroprosthetics by optimizing their functionality for individual patients.
A mathematical model makes it possible to describe the relevant characteristics of the patient’s skeletal muscle and to accurately predict the force as a function of the stimulation parameters. Indeed, synthesis of stimulation sequences or control strategies to achieve movement can be efficiently computed and optimized using numerical models [8]. Therefore, it can contribute to enhancing the design and function of controls applied to FES. Over the years, a great variety of muscle models have been proposed, differing in their intended application, mathematical complexity, level of structure considered, and fidelity to the biological facts. Some of them have attempted to exhibit the microscopic or macroscopic functional behavior, for instance Huxley [9] and Hill [10]. The distribution-moment model [11] constitutes a bridge between the microscopic and macroscopic levels. It is a model for sarcomeres or whole muscle which has been extracted via a formal mathematical approximation from Huxley crossbridge models. Models integrating geometry of the tendon and other macroscopic consideration can be found in [12]. A study by Bestel and Sorine [13], based on both microscopic Huxley and macroscopic Hill type model, proposed an explanation of how the beating of cardiac muscle is achieved through a chemical control input. It integrated the calcium dynamics in muscle cells that stimulate the contractile element of the model. Starting with this concept, we adapted it for striated muscle [14]. We proposed a musculotendinous model which considered the muscular masses and viscous frictions in muscle-tendon complex. This model is represented by differential equations where the outputs are the muscle’s active stiffness and force. The model input represents the actual electrical signal as provided by the stimulator in FES.
1.2. In Vivo Identification for Subject-Specific Parameters
In actual FES system, the appropriate tuning is achieved empirically by intensively stimulating the patient’s muscle for each task. If this adjustment could be calculated in the simulation, and if we could find the best signal pattern using virtual skeletal muscle, such method would be very helpful for movement synthesis for spinal cord injured (SCI) patients. However, in order to perform the simulation, an accurate skeletal muscle model is required to reproduce a well-predicted force for each muscle corresponding to the patient-specific characteristics.
For any biological systems, identification is a difficult problem due to the fact that (i) measurements must be as noninvasive, particularly on humans, (ii) some entities cannot be directly measured, (iii) an experimental setup and protocol have to be designed and certified, (iv) intersubject variations can be large, and (v) the large nonlinearity and complexity of the models cause some optimization algorithms to fail. Thus, few papers address biomechanical parameters identification in FES context, and they used macroscopic model for global force production [15]. Consequently, we described an approach for coupling the model with in-vivo measurements, that is, using a multiscale muscle model in an estimation procedure in order to perform the identification of the parameters, hence, giving access to physiologically meaningful parameters of the muscle. Preliminary result was reported in [16].
The skeletal muscle dynamics are in particular highly nonlinear, and we need to identify many unknown physiological parameters if a multiscale model is applied. The main objective of this paper is to develop an experimental computational method to identify unknown internal parameters from the limited information. For such identification of nonlinear system, extended Kalman filter (EKF) has often been used for skeletal muscle [17] and cardiac muscle identification [18]. It can work when the model is not complicated and not highly nonlinear. In this paper, a sigma-point Kalman filter (SPKF) was applied to the in vivo experimental data to identify internal states in the nonlinear dynamics of multiscale skeletal muscle model. SPKF has higher accuracy and consistency for nonlinear estimation than EKF theoretically. The feasibility of both identifications is verified by comparison. The computational performance is discussed.
The outline of the paper is as follows. The next section presents the formulation of the skeletal muscle model controlled by FES. The following section is devoted to the experimental identification of the model for isometric contraction, including the identification protocol. The experimental measurement was performed in-vivo on rabbits. Then, we present detailed results of the parameter estimation, the comparison with EKF, and cross-validation which illustrates the pertinence of the identification. Finally, we present some discussions, conclusions, and perspectives in the last sections.
2. Skeletal Muscle Model
Muscle modeling is complex, in particular when the model is based on biomechanics and physiological realities. Most of the muscle models have been based on phenomenological models derived from Hill’s classic work [10] and well summarized by Zajac [12]. Hill macroscopic model is a standard muscle model for practical use. Recent work [19] performed the validation of Hill model during functionally relevant conditions. They concluded that model errors are large for different firing frequencies and largest at the low motor unit firing rates relevant to normal movement. They pointed out that the reason may come from the Hill model assumption to consider muscle activation, force-length, and force-velocity properties independently. It was suggested that more physiological coupling between activation and force-velocity properties can be demonstrated in microscopic crossbridge models incorporating a dependence between physiologically based activation and crossbridge attachment [19]. Thus, our approach is to provide a multiscale physiology-based model on the both micro- and macroscale fact to obtain meaningful internal parameters.
Our muscle model is composed of two elements of different nature: (i) an activation model which describes how an electrical stimulus generates an action potential (AP) and initiates the contraction and (ii) a mechanical model which describes the dynamics in force (Figure 1). For details of the muscle model, refer to the previously published article [14]. Here, a brief summary of the model and the information necessary for the identification are given.
2.1. Activation Model
The activation model describes the electrical activity of muscle which represents the excitation-contraction phenomena. In muscle physiology, it is known that two input elements dominate muscle contraction: fiber recruitment rates and temporal activation [20]. The recruitment rates determine the percentage of recruited motor units. In Hill model, it has only one input which corresponds to this recruitment rate. Hill model is mainly applied for voluntary contraction where temporal activation is not dominant. However, in FES, since all muscle fibers are synchronously activated, this temporal activation which occurs by every stimulation pulse, is important. The recruitment rates can be determined on the pulse width (PW) and pulse amplitude of signal , generated by the stimulator [21]. The recruitment rate can be assumed to be a constant value when PW and remain constant.
The temporal activation can be considered as the underlying physiological processes which describes the chemical input signal, , that brings muscle cells into contraction as shown in Figure 2. Muscle contraction is initiated by an AP along the muscle fiber membrane, which goes deep into the cell through the T-tubules. It causes calcium releases that induce the contraction process when its concentration rises above a threshold and is sustained till the concentration drops back below this threshold once again. Hatze [22] gives an example of calcium dynamics () modeling. Since we focus on the mechanical response in this paper, we choose a simple model that renders the main characteristics of the dynamics. The contraction-relaxation cycle is triggered by the () associated with two phases: (i) contraction and, (ii) relaxation as in Figure 2. We use a delayed () model to take into account the propagation time of the AP and an average time delay due to the calcium dynamics. The frequency of the chemical input, , can be defined from the stimulation frequency. The time delay and the contraction time can be obtained from a twitch test by single stimulation pulse. A contraction takes place with a kinetics then, if no other AP has been received in the mean time, an active relaxation follows indefinitely with a kinetics . is linked to the rate of actine-myosine cycle whereas is related to the rate of crossbridge breakage. Finally, can be written as follows, where is a trapezoidal switching function which connects relaxation and contraction state from 0 to 1:
2.2. Mechanical Model
The model is based on the macroscopic Hill-Maxwell type model and the microscopic description of Huxley [9]. For crossbridge modeling, Huxley [9] proposed an explanation of the crossbridge interaction in a sarcomere. A sarcomere model can be used to represent a whole muscle which is assumed to be a homogeneous assembly of identical sarcomeres. The distribution-moment model of Zahalak [11] is a model for sarcomeres or whole muscle which is extracted via a formal mathematical approximation from Huxley crossbridge models. This model constitutes a bridge between the microscopic and macroscopic levels. Based on Huxley and Hill-type models, Bestel and Sorine [13] proposed an explanation of how the beating of cardiac muscle may be performed, through a chemical control input, connected to the calcium dynamics in muscle cell, that stimulates the contractile element of the model.
The model is composed of macroscopic passive elements and a contractile element controlled by input commands: the chemical input, , as suggested by Bestel and Sorine [13] for the cardiac muscle on the sarcomere scale and the recruitment rate, , on the fiber scale as shown in Figure 3. In order to express isometric contractions, whereas the skeleton is not actuated, our muscle model is introduced with masses (kg) and linear viscous dampers , () to ensure energy dissipation. On both sides of , there are elastic springs, , (), and viscous dampers to express the viscoelasticity of the muscle-tendon complex. The parallel element, , mainly represents surrounding tissues, but it can be omitted in isometric contraction mode. We assume the symmetric form with the two masses and passive elements identical. Here, , , and . and (m) are the lengths of and in the rest state. Initially, we can define the relative length variation as positive when the length increases, as in (2). In particular, in the case of isometric contraction, the following relationship exists (3): The dynamical equation of one of the masses is given by (4). and express the force and the stiffness of , respectively. The force at the output of this muscle model is the sum of the spring force and the damping force . When we measure the tension of skeletal muscle under in vivo conditions, the experimentally measured force is equivalent to . When we take the ratio of and , is offset and can be written as in (6). (7) is the relational equation in Laplace transform. From this relationship, the differential equation (8) can be obtained: Finally, in isometric contraction, the differential equations of this model can be described as follows: The dynamics of the contractile element correspond to (9) and (10). The terms and are the maximum values for and , respectively. From (3) and (4), the differential equation for is obtained as in (12). The internal state vector of this system should be set as .
3. Experimental Identification
In this study, we have developed a method to identify the parameters in the mechanical part of a skeletal muscle model. The input controls of the model are set as the static recruitment rate and the chemical control input from the activation model. These two controls are computed from the FES input signal. Note that the identification was performed with constant FES parameters for pulse width and intensity of electrical stimulation so that the recruitment rate is constant. The amplitude and pulse width were selected to recruit 90% of the maximum muscular force; then can be set as 0.9 for the fiber recruitment. In addition, the calcium dynamics in our model induce a time delay and an on/off control which represents contraction/relaxation so that correct data processing can avoid the detailed calcium dynamics. The trigger for signal can be calculated from the timing of the electrical stimulation. and are set as 20 and 15, respectively, as in [23].
In isometric contraction, the differential equations of skeletal muscle dynamics are straightly given in (9)–(12). In this case, , , , and are unknown time-varying values and , , and are unknown static parameters to be estimated. For the identification of this model, it is a nonlinear state space model, and many state variables are not measurable. Moreover, in-vivo experimental data include some noise. Hence, we need an efficient recursive filter that estimates the state of a dynamic system from a series of noisy measurements.
3.1. Sigma-Point Kalman Filter
For this kind of nonlinear identification, extended Kalman filter (EKF) is the well-known standard method. In EKF, the nonlinear equation should be linearized to first order with partial derivatives (Jacobian matrix) around a mean of the state. The optimal Kalman filtering is then applied to the linearized system. When the model is highly nonlinear, EKF may give particularly poor performance and diverge easily. In skeletal muscle dynamics, its state-space is dramatically changed between the contraction and relaxation phases. At this time, partial derivatives would be incorrect due to the discontinuity. Therefore, we introduced the sigma-point Kalman filter (SPKF). The initial idea was proposed by Julier and Uhlmann [24] and has been well described by Merwe and Wan [25]. SPKF uses a deterministic sampling technique known as the unscented transform to pick a minimal set of sample points (called sigma points) around the mean. These sigma points are propagated through the true nonlinearity. This approach results in approximations that are accurate to at least second order in a Taylor series expansion. In contrast, EKF results in only first-order accuracy.
An outline of the SPKF algorithm is described. For details, the reader should refer to [25, 26]. The general Kalman framework involves estimation of the state of a discrete-time, nonlinear dynamic system, where represents the internal state of the system to be estimated and is the only observed signal. The process noise drives the dynamic system, and the observation noise is given by . The filter starts by augmenting the state vector to dimensions, where is the sum of dimensions in the original state, model noise, and measurement noise. The corresponding covariance matrix is similarly augmented to an by matrix. In this form, the augmented state vector and covariance matrix can be defined as where is the state covariance, is the process noise covariance, and is the observation noise covariance.
In the process update, the sigma points are computed based on a square root decomposition of the prior covariance as in (15), where and is found using . is chosen in which determines the spread of the sigma points around the prior mean and is usually chosen equal to 0. The augmented sigma-point matrix is formed by the concatenation of the state sigma-point matrix, the process noise sigma-point matrix, and the measurement noise sigma-point matrix, such that . The sigma-point weights to be used for mean and covariance estimates are defined as in (16). The optimal value of 2 is usually assigned to : And these sigma points are propagated through the nonlinear function. Predicted mean and covariance are computed as in (18) and (19) and predicted observation is calculated as follows: The predictions are then updated with new measurements by first calculating the measurement covariance and state-measurement cross-correlation matrices, which are then used to determine the Kalman gain. Finally, the updated estimate and covariance are determined from this Kalman gain as below: These process updates and measurement updates should be recursively calculated in until the end point of the measurement.
3.2. Experimental Measurement for Identification
Stimulation experiments were performed on two New Zealand white rabbits at Aarhus University Hospital in Aalborg, Denmark, as depicted in Figure 4. Anesthesia was induced and maintained with periodic intramuscular doses of a cocktail of 0.15 mg/kg Midazolam (DormicumR, Alpharma A/S), 0.03 mg/kg Fetanyl, and 1 mg/kg Fluranison (combined in HypnormR, Janssen Pharmaceutica) [14].
The left leg of the rabbit was anchored at the knee and ankle joints to a fixed mechanical frame using bone pins placed through the distal epiphyses of the femur and tibia. A tendon of the medial gastrocnemius (MG) muscle was attached to the arm of a motorized lever system (Dual mode system 310B Aurora Scientific Inc.). The position and force of the lever arm were recorded. An initial muscle-tendon length was established by flexing the ankle to 90°. A bipolar cuff electrode was implanted around the sciatic nerve, allowing the MG muscle to be stimulated. Data acquisition was performed at a 4.8 kHz sampling rate. The muscle force against the electrical stimulation was measured under isometric conditions.
4. Result of Identification
In order to facilitate the convergence of the identification, the estimation process has been split into two steps. In the first step, we only estimate geometrical parameter . In the second step, we estimate the dynamic parameters: and . The biomechanical muscle model to be identified is presented as (9)–(12). We define the state vectors for geometric and dynamic estimations in SPKF as follows:
4.1. Parameter Estimation
Parameter estimation was performed using two rabbit experimental data. The stimulation signal input used for the estimation is composed of two successive pulses (doublet) at 16 Hz and 20 Hz with amplitude 105 μA and pulse width 300 μs. The stiffness has been estimated separately, using the experiments performed on the isolated muscle. The isolated muscle was pulled with the motorized lever. The stiffness is taken as equal to the slope of the straight line of the passive length versus force relationship in the experimental measurement as shown in Figure 5. The obtained was 4500 N/m. and can be obtained knowing the force response of muscle to a stimulation pattern with maximum signal parameter values (frequency, amplitude, pulse, and width). In this case, N/m and N were obtained from the measurement.
The experimental muscle force against the doublet stimulation (20 Hz) was used for parameter estimation during measurement updates for in SPKF. The covariance matrix has been initialized as in Table 1 for both geometric and dynamic parameter identifications. The evolutions of estimated internal states for and are obtained as shown in Figure 6. From these behaviors, we can confirm that the contractile element of the model is successfully tracking the contraction represented by the dynamics of differential equations under the estimation process. Figure 7 shows the estimated parameter and the error covariance. The fact that the error covariace is being decreased during identification process shows stable convergence of the estimation. Even with the different initial values setting, parameter estimation results converged into a particular value in SPKF as shown in Figure 7(a). We tested the estimation from 7 different initial values with an interval of 5 mm from 55 mm to 85 mm. The color is changed from magenta to cyan in the plot.
(a)
(b)
(a)
(b)
As can be seen from the resultant computational behavior on the graphs, the internal state vectors of the muscle model converged well to stationary values even with different initial values. After the complete estimation process for geometric and dynamic parameters, the estimated values are summarized in Table 2. Normalized RMS errors between measured and estimated force were also computed for each case. Moreover, the estimated length of the contractile element with each stimulation showed values close to the actual measured lengths of the extracted muscle. The border between the contractile element and tendon is not so clear visually, but it was approximately 7 cm for both rabbit GM muscles. The sizes of the two rabbits were similar; in particular, the estimated length of had good correspondence with the measured length for both rabbits and with different frequency stimulations. Physical parameters are impossible to make visual verification, but the estimated intrasubject property is maintained among the parameters obtained by different frequency stimulations from Table 2.
4.2. Comparison with the Extended Kalman Filter
The extended Kalman filter is generally chosen for nonlinear system identification. However, in EKF, first-order partial derivatives are used for the computation, which means that a matrix of partial derivatives (Jacobian) is computed around the estimate for each step. The detail of EKF is summarized in the appendix. When the process and measurement functions and are highly nonlinear, EKF can give particularly poor performance and diverge easily [25]. The estimate can have a bias due to the linear approximation especially for discontinuous systems.
Using the same computational conditions, such as initial values for states, parameters, and covariance matrices, parameter estimation was also performed with EKF to compare the estimation quality for this nonlinear system. The estimation results with EKF are shown in Figures 8 and 9. The numerical comparison for SPKF and EKF in 7 different initial conditions is summarized in Table 3 both for rabbit1 and rabbit2. With the observation noise covariance , the same as with SPKF, the result with EKF converged computationally as shown at the error covariance of EKF1 in the table. However, the matched result in force level does not mean directly that the estimated internal state and parameter are correct. The internal state in EKF estimation was not well estimated as the contracted strain variation. The difference can be observed by comparing the plots of Figures 9 and 6 for the EKF and SPKF estimates, respectively. In Figure 6, the estimated state of has two sharp peaks which correspond to the strain of the contractile element induced by muscle contraction by doublet stimulation. Thus, this internal state reflects the expected phenomenon in muscle dynamics. However, the estimated state in Figure 9 does not show sharp deformation induced by stimulation pulse.
(a)
(b)
The linear approximation in the EKF transformation matrix could be the cause of the lower quality of the state estimation. Finally, it resulted in a bias for parameter estimation. Thus, the converged value for under these conditions is not realistic. The EKF estimation was also performed with the observation noise covariance , giving more uncertainty to the measurement update. In this case, the estimated value became more realistic as the SPKF estimate did. However, the convergence became poor as in Figure 8(b) and it was highly dependent on the initial values as in Figure 8(a).
4.3. Model Cross-Validation
A cross-validation of the identified model was carried out to confirm the validity of this method. The resultant muscle force was calculated using the identified models with the control input generated by the stimulation frequency. Figure 10 shows the measured force response of the MG muscle of the rabbit and the simulated force with the identified model. The stimulation was three successive pulses in μA, μs, and Freq = 31.25 Hz. The red line indicates the measured muscle force, the blue and green dotted lines are the plots by the identified model with SPKF and EKF1, respectively. The model identified by SPKF could predict the nonlinear force properties of stimulated muscle quite well in this cross-validation. Figure 11 is the result in the case of stimulation trains with six pulses in 20 Hz. There is a good agreement between the measured and the predicted force; however, when muscle fatigue appears, the experimental force is lower than that from the simulation. Particualrly, the difference between the measured and the estimated forces at the end can be considered as the error coming from modeling without considering muscle fatigue. Apart from the error coming from muscle fatigue, we could confirm the feasibility of the FES muscle modeling and the effectiveness of the identification. Normally, the muscle property varies greatly, being dependent on the type of muscle, and the force response is highly subject-specific. Therefore, the muscle force prediction by cross-validation is quite difficult even for the approximate prediction especially when force production is predicted for each stimulation pulse. The identification based on experimental response would contribute strongly to realistic force prediction in electrical stimulation and FES controller development. Further investigation is required for the expression of fatigue phenomena [27], for its reproducibility and identification.
5. Discussion
The skeletal muscle model used in this identification protocol is based on both macroscopic and microscopic physiology which is unlike the black-box or other approaches using simple Hill muscle model. The structured model requires more parameters in highly nonlinear dynamics, which have to be estimated by experiment. However, the great advantage of our model is the insight which it can give to a muscle’s biomechanical and physiological connections, where the parameters have a physical significance such as length and mass. This paper describes an identification method which uses experimental response and a nonlinear, physiological muscle model to obtain subject-specific parameters. The advantage of animal experiment is to have direct access to extracted muscle for confirmation of the obtained parameter after the experiment.
This work was performed for application of FES stimulated muscles; however, since many living systems have nonlinear dynamics and subject specificity, this kind of identification approach itself could be applied to other organ models and clinical situations. In this work, both SPKF and EKF algorithms were applied for muscle dynamics identification. EKF showed a high dependency on initial settings of the computation and it was difficult to find the effective range of initial states and covariances. In SPKF, the obtained result was more consistent and robust with respect to various initial conditions. Moreover, for SPKF, the computation of a Jacobian is not necessary, so it can easily be applied even to complex dynamics. SPKF results in approximations that are accurate to at least second order in Taylor series expansion. In contrast, EKF results in first-order accuracy. Further, the identification accuracy is clealy improved, especially for nonlinear systems, but the computation cost still remains the same as for EKF. Advanced and robust system identification, including designing experimental protocols, has a very important role to improve the control issues in neuroprosthetics. In addition, the main difficulty in understanding human systems is caused by their time-varying properties. The systems are not static and change over time. The function of a human being is not always the same; for example, muscle fatigue can easily change the expected force response. In order to deal with the time-varying characteristics of a human system, robust biosignal processing and model-based control which corresponded to nonlinearity and time variance would provide a breakthrough in the development of neuroprosthetics [28].
6. Conclusion
We have developed an experimental identification method for subject-specific biomechanical parameters of a skeletal muscle model which can be employed to predict the nonlinear force properties of stimulated muscle. The mathematical muscle model accounts for the multiscale major effects occurring during electrical stimulation. Thus, the identification method was required to deal with the high nonlinearlity and discontinous states between muscle contraction and relaxation phase. The identified model was evaluated by comparison with experimental measurements in cross-validation. There was a good agreement between measured and simulated muscle force outputs. The result showed the performance which can contribute to the prediction of the nonlinear force of stimulated muscle under FES. The feasibility of the identification could be demonstrated by comparison between the estimated parameter and the measured value. In this study, the identification was performed by sigma-point Kalman filter and extended Kalman filter. The performance was compared and summarized under the same computational conditions. SPKF gives more stable performance than EKF. The internal state in EKF estimation was not well estimated as the contracted muscle strain. In SPKF estimation, keeping the realistic state transition and independence from initial conditions, it could realize converged solution for each identification trial.
SPKF is a Bayesian estimation algorithm which recursively updates the posterior density of the system state as new observations arrive online. This framework can allow us to calculate any optimal estimate of the state using newly arriving information. We believe that the proposed identification method has also the advantage for human muscle identification while it provides not only better accuracy in nonlinear dynamics but also adaptability to time-varying systems. The preliminary result is reported as in [29].
Appendix
Extended Kalman Filter
The extended Kalman filter [30] extends the scope of Kalman filter to nonlinear optimal filtering problems by forming a Gaussian approximation of the distribution of state and measurements using a Taylor series-based transformation. It is based on linear approximations to the transformation matrix. Assuming the same nonlinear system as in (13), EKF approximates a process with nonlinear difference and measurement relationships as follows: Note that the difference between EKF and KF is that the matrices and in KF are replaced with Jacobian matrices which are partial derivatives around the estimate in EKF. In order to avoid numerical problems, it is necessary to proceed the algorithm by using a square root factorization as a UD decomposition which ensures that the covariance matrix remains a positive definite matrix.