**Musculoskeletal Models in a Clinical Perspective**

Editor **Carlo Albino Frigo**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Carlo Albino Frigo Department of Electronics, Information and Bioengineering Politecnico di Milano Milan Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: www.mdpi.com/journal/applsci/special issues/ Musculoskeletal Models).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1818-3 (Hbk) ISBN 978-3-0365-1817-6 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

#### **Carlo Albino Frigo**

Carlo Albino Frigo, Department of Electronics, Information and Bioengineering, Politecnico di Milano, has been working in the area of movement analysis and motor rehabilitation since his graduation in 1976. His interests include human biomechanics, motor control and musculoskeletal models. He has been founder and president of the Italian Society of Clinical Movement Analysis (SIAMOC), the vice president of ISPO-Italy (International Society for Prosthetics and Orthotics), a member of the ESMAC (European Society for Movement Analysis in Adults and Children) committee, and a member of several scientific societies and editorial boards, participating to several European projects and educational courses. In 2012, he was an invited speaker at the Honorary Baumann Lecture at the ESMAC-2012 congress.

## **Preface to "Musculoskeletal Models in a Clinical Perspective"**

Musculoskeletal models are a promising tool for addressing many clinical problems in the next future. The increased computational power, availability of bio-images, and efficiency of the processing algorithms are all factors that make the musculoskeletal models more and more realistic and reliable. The scope of this book is to support the idea that musculoskeletal models can provide information not obtainable by other means and, thus, have an important role in the clinical processes. Bioengineers, medical physiatrists, orthopedic surgeons, physiotherapists, and movement scientists may be interested in the contents of this book and hopefully will be inspired to adopt the musculoskeletal modelling approach in their research area and to thus contribute to the evolution of this discipline.

> **Carlo Albino Frigo** *Editor*

### *Editorial* **Special Issue: Musculoskeletal Models in a Clinical Perspective**

**Carlo Albino Frigo 1,2**



**Citation:** Frigo, C.A. Special Issue: Musculoskeletal Models in a Clinical Perspective. *Appl. Sci.* **2021**, *11*, 6250. https://doi.org/10.3390/ app11146250

Received: 24 June 2021 Accepted: 30 June 2021 Published: 6 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

After the pioneering work of Scott Delp and colleagues dated 1990 (*An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures*, [1]) a huge amount of work aimed at creating a virtual representation of the human movement has been completed. The increased computational power and graphics capacity of computers have been fundamental in this process. Early in 2008, a very ambitious project was launched at the European level, with the objective of modelling all of the physiological phenomena of a human being: *The Virtual Physiological Human* [2]. The evolution of musculoskeletal modelling can be seen as a component of that ambitious project. Musculoskeletal modelling consists of reproducing the anatomical structure and the dynamic behaviour of a human subject, or a part of it, and simulating the effects of movement and forces upon the internal components. The interest in these techniques is very wide, spanning from sports biomechanics to ergonomics and clinical applications. In these different areas, the objectives and the methods can be different, but in all cases, musculoskeletal models can provide information that could not be obtained in other ways, considering, for example, the internal forces in a joint, the tension of the ligaments and the effects of muscle contractions. This Special Issue was particularly aimed at collecting experience and points of view concerning the potential of the musculoskeletal models in clinical applications. The review paper of Killen et al. [3] makes a point on the impact that musculoskeletal models and dynamic simulation have on clinical practice and in rehabilitation in particular. The perspective paper of BJ Fregly [4] provides a vision of what the role of musculoskeletal modelling in clinics could be and the steps required to make these techniques widely available and useful. Other papers report on interesting examples of problems in which musculoskeletal models have been fundamental to answer clinical questions: the complex mechanism of swallowing [5], where the role of different muscles is very difficult to investigate in other ways; the problem of lower back pain in pregnancy [6]; the question of the effectiveness of a particular reinforcement treatment for the abductor muscles in subjects affected by hip dysplasia after total hip joint replacement [7]. One paper [8] then provides an example of how the dynamic simulation can help us to understand the role of muscles during walking and answer a typical 'what if?' question: what happens if the Rectus Femoris muscle increases or decreases its activity in a specific phase of the gait cycle? Two more papers lay in the original mainstream of musculoskeletal modelling applications, that is, clinical orthopaedics. One faces the technical problem of understanding the effect of changing the tibial insert thickness in the knee joint arthroplasty [9], and the other reports on the forces occurring on the knee joint ligaments and bone surfaces while walking [10]. The papers collected in this Special Issue are obviously far from covering all of the possible applications of musculoskeletal models, but they do offer a demonstration of the usefulness of the modelling approach. This is truly a fast-growing discipline that is producing an enhancement of knowledge in many different areas. Being conscious of the potential effects of this evolution, which could be in terms of the prevention of musculoskeletal damage, treatment planning, or load assessment, is fundamental for improving the quality of health care services offered to people affected by neuromusculoskeletal diseases. This Special

Issue represents a contribution to the promotion of the continuous improvement and usage of this advanced methodology.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


*Review*

## **In Silico-Enhanced Treatment and Rehabilitation Planning for Patients with Musculoskeletal Disorders: Can Musculoskeletal Modelling and Dynamic Simulations Really Impact Current Clinical Practice?**

#### **Bryce A Killen \* , Antoine Falisse , Friedl De Groote and Ilse Jonkers**

Human Movement Biomechanics Research Group, Department of Movement Sciences, KU Leuven, 3001 Leuven, Belgium; antoine.falisse@kuleuven.be (A.F.); friedl.degroote@kuleuven.be (F.D.G.); ilse.jonkers@kuleuven.be (I.J.)

**\*** Correspondence: bryce.killen@kuleuven.be

Received: 31 August 2020; Accepted: 13 October 2020; Published: 16 October 2020

**Abstract:** Over the past decades, the use of computational physics-based models representative of the musculoskeletal (MSK) system has become increasingly popular in many fields of clinically driven research, locomotor rehabilitation in particular. These models have been applied to various functional impairments given their ability to estimate parameters which cannot be readily measured *in vivo* but are of interest to clinicians. The use of MSK modelling and simulations allows analysis of relevant MSK biomarkers such as muscle and joint contact loading at a number of different stages in the clinical treatment pathway in order to benefit patient functional outcome. Applications of these methods include optimisation of rehabilitation programs, patient stratification, disease characterisation, surgical pre-planning, and assistive device and exoskeleton design and optimisation. This review provides an overview of current approaches, the components of standard MSK models, applications, limitations, and assumptions of these modelling and simulation methods, and finally proposes a future direction.

**Keywords:** musculoskeletal modelling; musculoskeletal simulations; clinical research; rehabilitation; biomechanics

#### **1. Introduction**

The rise of high-powered computing has revolutionised both clinical and research fields, allowing the processing and collection of large amounts of data. In addition, new and complex methodologies have and continue to be developed addressing clinically relevant questions. These developments are present in many fields including cell biology, oncology, and drug development to name a few. One other field which has greatly benefitted is that of human movement biomechanics. The emergence of computational musculoskeletal (MSK) models in the past 30+ years has allowed great progress to be made in the field. MSK models present a mathematical description of the human musculoskeletal system (See Section 2 for detailed explanation), in terms of its segments, joints, and muscles as well as functional components including joint movement (i.e., kinematics) as well as muscle geometry and force-generating capacity.

The use of MSK models within both the research (considering all research fields, including sports, clinical, etc.) and clinical fields has and continues to rise (Figure 1). These physics-based modelling and simulation workflows provide an integrative platform that allows users to combine traditionally collected experimental data (i.e., three dimensional (3D) motion capture) with information on 3D structure and function of the MSK system (i.e., MSK model) allowing more advanced parameters to

be estimated. Specific parameters which can be estimated using these methods include joint angles (i.e., kinematics) and moments (i.e., kinetics), but extend to MSK structure loading (i.e., ligament and muscle forces) and joint contact loading. Many of which cannot be estimated without the use of MSK models. In addition to studying these parameters in isolation, interactions (e.g., effect of kinematics/moments on joint contact loading) can also be investigated. These additional parameters can be applied to various MSK conditions and numerous stages of the clinical treatment pathway.

**Figure 1.** Total numbers of papers (blue) using musculoskeletal (MSK) modelling (Source: PubMed https://pubmed.ncbi.nlm.nih.gov/), and number of clinical science papers (orange) using MSK modelling (Source: Dimensions https://app.dimensions.ai/).

The remainder of this paper will present three sections: first, an educational component explaining MSK models (See Section 2.1) and their associated simulation methods (See Sections 2.2–2.8). The second section provides examples of specific clinical applications where MSK models have been and are currently being introduced to enhance treatment. Specifically, using MSK modelling and simulations, MSK diseases and conditions can be investigated to determine disease-specific characteristics (see Section 3). Furthermore, by testing patients at multiple time points, disease-specific biomarkers can be determined which predispose patients to disease progression, giving insights into disease aetiology (see Section 3). In the context of MSK rehabilitation using MSK modelling and simulation tools, the precise load applied to MSK tissues during functional activities can be estimated and compared to loading guidelines. Physical therapy or exercise regimes can then be modified to provide optimal, patient-specific MSK loading to deliver the optimal regenerative response (see Section 4). In addition to the investigation and optimisation of MSK loading, MSK modelling and simulations can be employed in the design and optimisation of active prosthetics to improve motor function. Likewise, critical parameters for pre-surgical planning procedures, implant and assistive device design can be evaluated and their impact on post-intervention functionality evaluated *in silico* (see Section 5). The applications introduced in this review have a strong emphasis on lower-limb joint contact loading due to the high prevalence in literature and its widespread applications to various MSK conditions. The third section of this review presents a discussion on modelling and simulation considerations and a roadmap for the future integration of MSK modelling and simulations in clinical practice. Indeed, despite clear integration pathways and advantages for clinical applications, the use of

MSK modelling and simulations remains scarce in standard clinical practice. Poor integration is likely multifactorial with barriers, including lack of easy point-of-care integration, specialist knowledge for use and interpretation and, in some cases, lack of knowledge of the added benefits of MSK modelling and simulations (see Section 8). Additionally, technical limitations exist regarding the accuracy of simulation results, simplifications made within MSK models, and the often cost- and time-intensive frameworks (see Section 7).

#### **2. Model-Based Analysis of MSK Function**

#### *2.1. Model Structure*

Musculoskeletal models are, at their core, a mathematical representation of the MSK system (Figure 2). Here, we will describe the general structure of commonly used MSK models. On the simplest level, models contain simplified representations of bones, connected by a series of joints which are then actuated by muscle tendon units (MTU) typically represented as line actuators that have force-producing properties.

**Figure 2.** Left: Example of generic full-body MSK model [1] detailing each model component. Displayed are each of the bones, an example of an experimental marker set (pink dots), and the biceps femoris long head, vastus medialis and medial gastrocnemius muscle tendon units. Right: Example of a more complex model [2,3] with representations of ligaments and cartilage surfaces of the tibiofemoral and patellofemoral joint.

Each bone of the skeletal system is represented as a rigid body with a geometry, a mass, and inertial properties.

Joints connect two adjacent bodies and define their articulation (i.e., location and orientation of joint axis around which the body moves) and permitted motions (i.e., degrees of freedom (DOF)). Several different joint definitions are used within MSK models which differ in the number of DOFs (i.e., movement directions). These range in complexity from allowing no movement (weld joint e.g., fibrous joint between the tibia and fibula), permitting a single rotation (hinge joint, e.g., dorsi-/plantar-flexion of the ankle), allowing three rotations (complex ball and socket joint, e.g., hip), all the way up to joints permitting three rotations and three translations (6 DOF joint, e.g., tibiofemoral). Joint kinematics describe how bony and soft tissues limit joint movement (e.g., defining hip rotation range of motion due to the hip capsule). As a result, the permitted range of motion of each DOF is defined, representing these bony and soft tissues structure function, without explicitly modelling

them. Some models introduce more complex joint definitions, (e.g., through kinematic coupling) which relates an undetermined (i.e., difficult or impossible to measure) DOF to a measured DOF. A specific example of this kinematic coupling is the patellofemoral joint [4]. Given the difficulties in reliably tracking its movement with markers, patella movement is defined as a function of the knee angle rather than based on measured marker positions.

Muscle tendon units span joints and actuate them by developing force along their line of action (i.e., path). In MSK models, muscles and tendons are modelled as a combined unit, hence the term muscle tendon unit. The path of each MTU is defined by an origin and insertion point, but also the intermediate pathway can be specified to represent wrapping over soft tissue or bony structures. This intermediate pathway can be defined using a combination of static or moving points, or analytical surfaces [5]. The pathway determines the elongation of the muscle as a function of the joint movement and by extension the moment arm [6]. The force-producing capacity of the MTU is modelled based on a Hill-type model of muscle contraction [7]. This model describes generic curves describing force generation of active and passive MTU components as a function of length and velocity, thereby capturing the effects of length and contraction velocity on force generation. These generic force-generating curves are made actuator specific by assigning MTU-specific parameters (e.g., maximal isometric force, optimal fibre length, tendon slack length, pennation angle and maximal fibre contraction velocity) to appropriately scale and tune the force-generating capacity.

These MSK models can then be used in combination with 3D experimental gait data (i.e., marker positions, and ground reaction forces) to estimate biomechanical parameters during dynamic function—for example, joint angles—but also muscle activations and forces ultimately allowing calculation of joint contact forces. Such a workflow is referred to as an inverse modelling workflow and consists of the following different steps (Figure 3).

#### *2.2. Model Customisation*

First, a customised MSK model needs to be generated. This is typically done by scaling a generic model to match the subject's gross anthropometry or by utilising subject-specific geometry obtained from medical imaging.

Scaling of a generic model is performed with the aim of matching the gross anthropometry of a subject. This process linearly scales each body (e.g., femur and pelvis) to match measured segment lengths, anthropometric table values, or most commonly, experimentally acquired motion capture marker positions [8–10]. Along with model geometry, other parameters including joint positions, segment mass and inertia, and MTU pathways are also scaled.

Alternatively, personalised models can be used [11–15]. Using medical imaging segmentation, bone and joint (i.e., cartilage and ligament) geometries and positions, and MTU pathways can also be personalised [11,16,17]. Different types of medical imaging allow the inclusion of different personalised parameters, for example, using computer tomography (CT) information regarding bone geometry can be obtained, however for information regarding soft tissues (i.e., muscle, ligaments, and cartilage) magnetic resonance imaging (MRI) is required. In the absence of complete medical imaging segmentations, statistical shape modelling methods can be utilised to reconstruct complete subject-specific geometry from incomplete data [16,18–20]—see Section 9. The use of 2D ultrasound (US) can provide information regarding MTU parameters (e.g., tendon resting length), and combined with motion capture (i.e., 3D US), information regarding MTU and bone geometry can be obtained [21–26]. The collection of muscle strength data, via dynamometry allows for the customisation of the maximal force-generating capacity (i.e., maximal isometric force) of muscle groups. Further, these dynamometry data can be used to further refine parameters within the previously mentioned Hill-type muscle model [26–29]. Despite the vast number of parameters which can be personalised, the inclusion of the above personalised parameters is uncommon and building fully subject-specific MSK models is rare given the complexity of the additional data collection. Typically, a selection of these personalised data is used to customise existing generic models.

#### *2.3. Inverse Kinematics*

Once the MSK model has been customised (either through scaling or personalised methods), 3D marker positions acquired within a gait laboratory based on a motion capture system can be used to estimate the joint angles, velocities and accelerations (i.e., joint kinematics) throughout the movement (e.g., walking, running, stair climbing) using an inverse kinematics procedure. In such procedures, the error between recorded marker positions and model marker positions is minimised by varying the segment positions in accordance with the permitted DOFs (i.e., rotations or translations) and imposed constraints [30,31]. The resulting analyses provide an estimated joint angle for each model DOF, at each time point of the specified task. From these kinematic analyses, additional parameters such as MTU length and moment arms can also be estimated. The force-generating

capacity of the MTU can be calculated, accounting for its length and lengthening velocity, as derived from joint kinematics. Combined with MTU moment arm information, the potential of the MTU to induce segmental acceleration and therefore its contribution to the movement can then be evaluated (i.e., the moment generating capacity of the MTU).

#### *2.4. Inverse Dynamics*

In a next step, model and estimated kinematics can then be combined with ground reaction forces acquired via force plates to estimate the joint moments during the same movement. To estimate joint moments at each joint, and each DOF, typically an inverse dynamics approach [32] is taken using a classic mechanical formulation. The joint moments represent the resulting action in the joint and are representative of the gross muscle action required during movement.

#### *2.5. Muscle Activation and Force Estimation*

Following inverse dynamics, individual muscle activations and forces can be estimated which balance the joint moments. Muscle activations are typically estimated using optimisation-based methods that use mathematical criteria to divide the required muscle forces at a joint between all muscles spanning that specific joint. The summed muscle forces must satisfy the constraints of the system (i.e., produce the required muscle moments to balance the joint moments subject to an objective function). Typical objective functions used within the literature are the minimisation of muscle activations to the *nth,* typically the second, power. Muscle activations and forces can be solved using various optimisation approaches, most commonly static optimization [33–35]. More complex implementations include dynamic optimisation [34,36–38] which can account for more complex muscle tendon dynamics. Another method is computed muscle control [39,40] which couples static optimisation with feedforward and feedback controls to approximate a set of desired kinematics (e.g., from inverse kinematics) and simultaneously estimate both muscle activations and forces. These optimisation methods can be supplemented with experimentally acquired electromyography (EMG) data to constrain estimated muscle activations [41,42]. Calculated muscle activation and force patterns allow researchers and clinicians to imply any changes in muscle function due to different pathologies, or locomotion strategies.

#### *2.6. Joint Contact Force Estimation*

With knowledge of muscle forces, analyses of joint contact forces can be performed. Typically, MSK models are rigid body models, meaning ligaments and cartilage surfaces are not considered. The output of a typical joint level analysis is therefore an intersegmental contact force [43] considering the joint reaction forces, muscle forces, and body velocities and accelerations. In contrast to finite element analysis (FEA) that provides local tissue loading and deformation, joint contact forces are representative of the overall joint loading given the effect of the forces through foot–ground contact, and the effect of muscle forces. While these approaches give an appreciation for the forces between two bodies, typically these are used as a surrogate for the true variables of interest, cartilage on cartilage contact forces/pressures that can be calculated using more complex models (see below) or different modelling techniques (i.e., FEA).

#### *2.7. Complex Models and Modelling Approaches*

More advanced models have been developed that include additional components, such as cartilage and ligament surfaces [3,44], allowing the estimation of ligament strains, and forces, along with cartilage on cartilage contact (Figure 3). Estimations of cartilage on cartilage contact has been performed, in MSK models, typically using either an elastic foundation [3] or Hertzian contact model [45], where force is estimated based on overlap or penetration of two surfaces. In addition to cartilage-level analysis, these more complex models allow for ligament level analysis. The inclusion of ligaments in MSK models allows for the estimation of ligament strain and subsequently forces from inverse kinematics

analysis (i.e., due to changes in joint angle). Furthermore, estimates of ligament forces can be used in the estimations of joint kinematics in complex six DOF joints. Such methods have been applied to the knee joint [3] where estimations of joint kinematics are subject not only to the marker positions, but also cartilage on cartilage contact and ligament forces. These complex MSK models, considering both cartilage and ligaments, provide models and simulation results which can better represent *in vivo* joint mechanics.

#### *2.8. Forward Simulations*

An alternate method to the inverse approach detailed above are forward simulation methods which use muscle activations and the resulting muscle force to drive MSK models and subsequently estimate the resulting kinematics and dynamics. As such these models allow users to establish a causal relationship between the forces acting on the human body and the resulting movement, thereby allowing an analysis of 'what if' scenarios. An extension of these forward simulations are predictive simulations, which can be employed to generate novel movements without relying on experimental data. Predictive simulations are typically used to identify muscle excitations and therefore movement patterns that optimise a user-defined cost function (e.g., minimise joint contact forces, muscle fatigue, etc.) subject to constraints describing, among others, the dynamics of the MSK system. These methods exist stand-alone [46] but have also been integrated into existing software packages [47,48].

The following sections will detail specific examples of the use of MSK modelling and simulations and detail pathways for integration into clinical practice. A majority of the detailed applications are focused on lower-limb MSK loading due to the extensive available literature and its relevance to a range of MSK conditions rather than a specific pathology. However, these methods can and have been extending to other segments including the spine [49], upper limbs [50], and head–neck [51].

To conclude, parameters and analyses made available through the use of MSK models and modelling provide valuable information that can be used for not only research studies on clinical populations, but also enhancing clinical treatment pathways, as described below. Indeed, MSK models combined with the tools and methods described above allow detailed analyses of muscle and joint contact forces, in addition to joint angle and moment analysis that most clinical practice software (e.g., Vicon Plug-in-gait (Vicon Motion Systems, Oxford, UK)) provide. Facilitating the adoption of these models are a number of software packages, both commercial (e.g., AnyBody [52]) and open source (e.g., OpenSim [53,54]). Both modelling and simulation platforms are typically released with different generic models and are supported by a combination of dedicated developers and their respective research community. The provision of predictive simulation techniques was just recently facilitated through the extension of commonly used software (i.e., OpenSim) to include these predictive simulations as a standard tool [47,48].

#### **3. Role of Musculoskeletal Models in Disease Prevention and Patient Stratification**

The use of MSK modelling and simulations in cross sectional studies (i.e., single time point) can provide valuable information into disease states and biomarkers (e.g., joint contact forces) that cannot be (easily) measured *in vivo*. For example, MSK modelling and simulation methods have been applied to patients who have undergone knee ligament reconstructions [55] who are at an elevated risk of osteoarthritis (OA) development. This was previously thought to be driven by higher joint contact forces following reconstructive surgery, it was shown, utilising MSK modelling and simulations, that these patients exhibit lower joint contact forces. This suggests that lower knee joint contact forces, instead of higher, play an important role in OA development. Likewise, others [56] were unable to indicate significant differences in knee joint contact forces in subjects with focal cartilage defects, following traumatic injury. Based on this insight, the magnitude of joint contact forces may be used as a biomarker to monitor disease development in this specific population. These biomarkers can then inform targeted rehabilitation, be used as model-based outcome measures, and as early warning signs or detectors for developing diseases which do not present in typical clinical tests or static medical

imaging modalities. Likewise, the use of MSK modelling and simulations enables researchers to determine differences in the role of altered mechanical joint environment in OA at different lower-limb joints (i.e., knee vs. hip OA). For example, during gait, joint contact forces in the involved joint is elevated in knee OA [57], whereas elevated joint (i.e., hip) contact forces are not present in hip OA patients [58,59].

Although MSK modelling and simulations are highly valuable in investigating MSK diseases, a potentially more powerful use is in preventing and slowing the progression of different MSK conditions. In this respect, MSK models can be employed within longitudinal trials to study changes in joint contact loading over time [60] (Figure 4). These longitudinal studies allow researchers/clinicians to identify patients who do and do not develop/progress with respect to a specific MSK condition. Examining parameters such as joint moments and in particular joint contact forces would allow researchers/clinicians to identify potentially modifiable parameters to slow disease progression. These biomarkers can then be targeted in rehabilitation such as gait retraining to reduce the likelihood of disease development and limit disease progression. While the above-mentioned longitudinal monitoring has been carried out in typically small cohorts, the type of analysis mentioned above (i.e., biomarker identification) is yet to become mainstream. One of the major limitations is the absence of longitudinal movement analysis data in large longitudinal studies. Including 3D motion capture data in the follow-up of OA patients, together with, medical imaging and various questionnaires which assess subjective disease state, would allow untangling the role of mechanical joint loading in disease progression. With increased computational speed and adequate data collection, these longitudinal studies will become more common in the future as simulation methodologies are more accepted by the medical community.

**Figure 4.** Longitudinal changes in first and second peak hip joint contact forces in patients before, and 3 and 12 months following total hip arthroscopy [60].

#### **4. Role of Musculoskeletal Models in Rehabilitation**

Along with investigating MSK disease state, and aetiology, MSK modelling and simulations can provide benefits relating to MSK rehabilitation. Rehabilitation is typically aimed at restoring muscle strength, joint stability, or joint loading to a typical or pre-injury state. Prescription of rehabilitation is typically based on general guidelines and are not highly tailored to individual patients. It is in this generic prescription where MSK modelling and simulations can enhance the rehabilitation pathway and decision making. By investigating MSK loading during various rehabilitation exercises and determining the load applied to specific tissues, exercise dosage can be tailored to achieve the optimal loading and inform rehabilitation progression [61,62]. Previous research [63] in this area has provided insights into the magnitudes and locations of knee joint loading during different rehabilitation exercises (Figure 5). This information can be particularly useful to personalise the prescribed exercises based on injury location and rehabilitation stage. Utilising the disease biomarkers mentioned previously, the effectiveness of different rehabilitation programs can also be assessed in terms of restoration of typical joint loading or effect on a targeted biomarker. In addition to the investigation and refinement of current clinical guidelines, MSK modelling and simulations can be employed to, for instance, evaluate the effectiveness of new rehabilitation exercises to restore appropriate loading. For example, based on MSK modelling and simulations input, a ranking of exercise regimes was proposed to prevent femoral neck bone loss in females at risk of osteoporosis [64,65]. Using MSK modelling and simulations, compensatory movement strategies to lower knee joint contact forces in patients with knee OA during stair climbing could be identified [57]. These insights may assist clinicians in providing guidelines for rehabilitation of patients to potentially avoid further progression and joint degeneration.

**Figure 5.** Average tibiofemoral joint contact forces for a range of different functional and rehabilitation tasks [63].

This is even more clear in the field of gait retraining, where the therapist instructs the patient to modify their walking pattern in an attempt to alter the loading on a specific MSK tissue, typically cartilage. These retraining methods have been initially applied to patients with OA whereby specific kinematic gait features (e.g., toe position, trunk lean) are modified to reduce the moments around the knee [66–71] or hip joint [72]. Similar to rehabilitation exercises, these gait retraining goals are based on generic guidelines and the effect on the MSK loading in the individual patient is typically not verified. Over the last decade, the optimisation of MSK modelling and simulation frameworks has driven the development of real-time feedback tools for gait retraining [73–76]. These tools allow clinicians to monitor gait patterns and estimate joint moments in real-time to assess their adherence to the prescribed gait strategy but also to evaluate its effect on MSK loading. Although the testing of different gait strategies can be performed based on experimentally acquired data of the patient's gait, newly developed predictive simulations have been employed to test the effect of an imposed gait strategy on MSK loading without the need for repetitive and time consuming in-person patient assessments. The effect on various MSK tissues can then be investigated to determine if the gait strategy is a viable rehabilitation option. If further developed, this may provide not only optimised treatment planning and likely improved rehabilitation, but also significant time and cost savings for both the patient and clinician. In conclusion, incorporation of MSK modelling and simulations in the evaluation of rehabilitation programs can provide a critical evaluation of the MSK effect of current and proposed rehabilitation strategies and can thus be used to determine if a long term effect of a dedicated program is warranted. Such investigatory studies are currently not standard practice, however the use of MSK modelling approaches to evaluate relevant parameters (i.e., joint contact loading) should be considered for future research.

#### **5. Role of Musculoskeletal Models in Pre-Surgical Planning, and Implant and Assistive Device Design**

Previously mentioned clinical treatment applications were primarily focused on physical therapy and rehabilitation, however MSK modelling and simulation applications also extend to orthopaedics. The testing of joint implants has long been performed in the engineering field to improve patient satisfaction, implant lifespan and design, and post-surgical functionality. Implant testing is now often performed using standardised loading procedures, but in the future the use of *in silico* approaches, performing virtual implant testing will become part of the standard design and testing process enabled and enhanced by MSK modelling and simulations. Using MSK modelling and simulations, population-based or even patient-specific loading parameters (i.e., magnitude and orientation) can be used to perform *in silico* implant testing, as well as expand on current experimental testing rig protocols. The inclusion of physiological loading data is highly relevant and the logical next step especially in implants which were designed using patient-specific data [77–83]. Another factor also related to the one-size-fits-all design of custom implants, MSK can help in providing guidance to the physical design or alignment of implant components (e.g., femoral head and acetabulum cup positioning—femoral neck off set). Previous studies [84], albeit in native knees, has illustrated the coupling between tibial alignment and knee joint contact force magnitude and location, as well as ligament strains. These insights, in addition to illustrating the importance of tibial component alignment, may serve as guidelines for component alignment during surgery, especially given the importance on soft tissue balancing. Similarly, these assessments have been performed in total knee joint replacements [85–87], investigating the relationship between component alignment and loading parameters (i.e., component stress). Furthermore, using MSK modelling and simulations, identification of gait kinematics which increase the risk of edge loading following total hip replacement could be determined, thereby identifying movement patterns to be corrected during the post-operative phase [72]. Using MSK modelling and simulation workflows, it is already feasible to study the interaction between patient gait characteristics, implant type, and surgical approach, allowing surgeons to optimally tune the surgical intervention based on function rather than patient geometrical features [88]. Ultimately, both design (i.e., sizing, component radius) and alignment (i.e., stem alignment and head angle) parameters can be assessed *in silico* to design patient-specific implants based on the previously described predictive simulations. This would allow fully optimised post-operative functionality not only at the level of implant loading, but also at the level of gait kinematics and maybe, one day, energy expenditure.

Similar design optimisation (i.e., arch size, location, and manufacturing material) can be applied to shoe insoles to study the effect of variables such as joint contact forces and plantar pressures [89–93] and, more recently, in the optimised design of insoles [94]. This approach is highly relevant in the area of diabetic foot care, where individual insole design is key in preventing foot ulcer formation [89,95]. Further, MSK modelling and simulations can be used in the optimisation of assistive devices (Figure 6) and exoskeletons. In addition to the design parameters discussed previously, control of active components (i.e., motor torques) in these devices can be optimised using MSK models and simulations [94,96–99]. With respect to exoskeletons, the use of predictive simulations can be used to optimise control strategy (e.g., for balance control or evaluation of comfort based on contact pressure evaluation [100,101]).

**Figure 6.** Example of a MSK model implemented in OpenSim which includes a model of a transtibial prosthetic device [102].

#### **6. Model Validity**

As MSK models and their associated modelling and simulation methods estimate parameters which cannot easily be measured *in vivo*, direct validation is not always straightforward. As such, various simulation standards are set out within the MSK modelling and simulation community to ensure feasible and valid solutions. These include assessment of errors between acquired and modelled marker trajectories, and assessment of residual forces following inverse dynamics [103]. More advanced validation steps include comparing estimated muscle activations to acquired EMG data to illustrate muscle on–off timings and muscle activation magnitudes. These comparisons are done using correlations where current methods show correlations ranging from 0.46–0.99 [104]. Along with these stand-alone validation steps, some initiatives have provided validation data sets for joint contact forces and joint kinematics derived from instrumented endoprosthesis [105–108]. Using these validation datasets, estimation of joint contact loading magnitude within 0.1–0.5 [3,109] and 0.5–1.0 [110,111] body weights for the knee and hip joints, respectively, have been reported. These studies, although showing good agreement with validation data, cannot be carried out for each subject in large cohort studies. A number of studies using these datasets have already provided not only confidence in the modelling and simulation approaches but important methodological advances and consideration [109,112,113]. With such datasets available, and becoming more common, the refinement of methods, and the validity of simulations will continue to improve.

#### **7. Modelling and Simulation Assumptions**

Although MSK models provide the opportunity to investigate various pathologies and *in vivo* parameters, just as important as their strengths, limitations are present and must be considered. These limitations relate to both modelling and simulation assumptions and implementations.

A commonly criticised assumption within MSK models is their simplified representation of MTU parameters, and the activation/force estimations. This stems from two origins, first the generic MTU parameters being estimated from cadavers. To partially overcome this, numerous calibration and optimisation techniques have been proposed [41,113–118] and implemented to produce a set of MTU parameters which are likely more accurate than those in generic estimates. Second, the previously mentioned optimisation methods, although providing a set of forces and activations which balance the MSK system, may not accurately represent various features including antagonistic muscle contraction, and muscle spasticity in specific pathologies such as cerebral palsy. With respect to muscle activations, estimates can be improved by using experimentally acquired EMG [113,119]. Additionally, further processing and pathology-specific procedures have been developed to better represent known differences in muscle function [118]. The sensitivity of simulations to differences in MTU parameters and optimization method have shown variations in estimated hip joint contact force from 2.5 body weights to 5.4 body weights [111].

Along with these MTU related assumptions, simplified representations of anatomy may affect simulation result accuracy. For the lower limbs, the ankle–foot is often represented with a single plantar-/dorsi-flexion DOF despite the numerous joints and DOF present *in vivo*. This simplification has implications when calculating model parameters [120]. Additionally, the tibiofemoral joint is often represented as a single DOF hinge joint where kinematic constraints are often used for the remaining DOFs. To overcome these joint model simplifications a number of methods have been proposed to more accurately represent *in vivo* joint motions [121–126]. These simplifications also exist in the upper body where the torso is typically represented as a single unit without muscles and highly simplified articulation not representative for the complexity of spinal movement [49].

#### **8. Barriers for Clinical Implementation**

The adoption of MSK models, modelling, and simulations into clinical research and practice is hindered by additional barriers, independent of the above modelling and simulation assumptions. The first barrier is knowledge of the benefits of MSK modelling and simulations within clinical practice and the added value and information they can provide. Although within research fields and journals advancements are well documented, transfer of knowledge to relevant clinical partners is often scarce. As such, clinicians may be unaware of cases where MSK models can provide fast and reliable answers to clinically relevant questions. Therefore, there is a clear need for education of the various applications.

Along with education, simple implementation of complicated workflows should be developed and made available to clinicians. These workflows should be able to provide the desired parameters (e.g., joint moments, joint contact forces) without the need to execute complete time-consuming workflows. Such approaches [127,128] have been implemented for simple applications where gastrocnemii and gastrocnemius muscle lengths are estimated using MSK modelling and simulation outputs as reference and clinically feasible measures (i.e., kinematics). Additionally, specific clinically orientated software has been developed for the analysis of gait [129] where simplified approaches have been used, with previous research showing similar accuracy to more complex research driven software [130]. However, the extension of similar interfaces to estimate more complex parameters remain largely unavailable, although some progress is being made in the area of pre-operative planning in cerebral palsy [131]. Along with these easy to use applications, training for both data collection and simulation workflows is imperative for clinical implementation. Additionally, training and simplification of workflow output, will ensure simulation results can be implemented in clinical practice with minimal burden. To ensure MSK simulations are valid and reliable, standardisation of both models and methods is required. A number of different models and frameworks are used throughout research to address different questions; however, to allow comparison of simulation result both between patients and time points, standardisation of modelling frameworks is imperative.

#### **9. Roadmap for Future Transfer of MSK Modelling and Simulations into Clinics**

The final section of this review provides a roadmap for future applications, developments, and improvements to current MSK modelling and simulation workflows. Several application fields are quickly evolving and will likely provide a next generation of treatment pathways, which build on current applications. Specifically, surgical related applications, both in the design of prosthetics and surgical technique optimisation. Personalised representations of patients, coupled with different implant/prosthetics designs allow for virtual device testing (i.e., loading of the implant/prosthetics). These virtual testing can be performed using a combination of pre-surgical data and predictive simulations. Utilising these predictive simulations further, the effect of different surgical techniques (e.g., femoral de-rotation, tendon transfer) can be refined to achieve the optimal treatment outcomes. If these MSK modelling and simulation methods are to be extended to those requiring surgery, accuracy in modelling not only otherwise healthy individuals but those with pathologies is imperative. A specific challenge is appropriately modelling those with motor impairment, and MSK disorders. These may relate to altered muscle activation strategies, bony deformities, or altered MTU characteristics (e.g., muscle weakness). Many of these applications, particularly surgical planning, hinge on highly personalised representation of the MSK system, specifically through subject-specific MSK models.

Despite a number of researchers and groups developing these subject-specific models and modelling platforms, their use especially in clinical research is scarce. This is due to the large amount of data which is required to build these models, along with the cost and time burden to collect this data. Simplification of this personalisation process is required if these models, and applications are to be realised. With respect to anatomical data (i.e., bone segmentations), the use of incomplete or sparse data to synthesize full data is a promising pathway. To this end, the use of statistical shape models (representative of the average shape of a specific population) and associated principal components (describing the shape and size variation) is becoming more common. These methods can use either incomplete segmentations [18,19,132] or a motion capture trial [18,133,134] to reconstruct personalised lower-limb geometries and assemble them into a MSK model. While becoming more common, the application of similar methods to other MSK tissues (muscles, ligaments, and cartilages) would allow for the generation of fully personalised MSK models from low-fidelity clinically feasible data. An alternative approach here is the use of population-based modelling approaches to evaluate model predictions for a range of parameters rather than aiming to accurately represent the individual patient [135,136].

Along with the need to use low-fidelity data to create these models, low-fidelity data for executing these workflows is also required. For example, reducing the need for 3D motion capture, force plate, and EMG data will facilitate the use of MSK models in clinical practice. Using inertia measurement units (IMUs) combined with machine learning and artificial intelligence methods, joint angles, moments, and joint contact forces [137–143] can be estimated in standard clinical practice. These machine learning methods have also been applied to estimate gait speed, cadence, knee flexion angle, and common gait impairment metrics using 2D video [144]. Joint kinematics estimated using such methods could potentially serve as low-fidelity inputs into the previously discussed MSK modelling frameworks. These machine learning methods, combined with low-fidelity data (e.g., IMU, 2D video) will likely provide the most practical clinically feasible assessment of gait. Likewise, collection of a subset of muscle activity may be sufficient to provide an estimated of the muscle coordination strategy without the need of collecting a large number of EMG channel. These data can then be used to improve estimates of muscle forces and activations. Similar to the previously mentioned statistical shape modelling methods, using a subset of collected EMG, muscle activations can be reconstructed through the use of muscle synergy methods [145–149]. However, the use of these low-fidelity data methods, combined with simplified MSK modelling for simulation purposes still needs to be further explored in a research context prior to its adoption in clinics.

In this respect, the adoption of low-fidelity input data for both model creation and execution should however be carefully coupled with both uncertainty and sensitivity analyses to evaluate the set of parameters that have a large effect on simulation results, and as such should be accurately defined. Conversely parameters which have little effect on simulation results can be defined using low-fidelity and less accurate methods, hence reducing the data collection burden. Such analyses have already been performed showing the influence of mass and inertial properties [150]. joint location estimations [123,150], bone–joint geometry and alignment [84,151], muscle parameters [38,152,153], and muscle redundancy formulation [111] on various simulation results including joint moments, muscle force estimates, and joint contact forces. Despite these previous studies, these forms of sensitivity and uncertainty analyses have not been systematically completed on a single dataset, meaning the ranking of these parameters, in terms of simulation sensitivity is difficult to interpret. Additionally, previous research has typically investigated the effect on a single or limited simulation results. As the

importance of a specific parameter is likely dependent on the output of interest, extensive sensitivity analysis on an output by output basis will provide the ideal balance between data collection burden and simulation validity.

In conclusion, the integration of clinically collected data (i.e., 3D motion capture) with MSK models can provide estimates of numerous clinically relevant parameter related to MSK loading which cannot be estimated using commercially available clinical gait analysis suites. Estimates of these biomarkers can enhance several clinical treatment pathways including rehabilitation, patient stratification, and disease aetiology investigation. Additionally, optimisation of endoprosthesis, prosthetics, assistive devices, and exoskeletons design and controls can be performed using novel predictive simulations. Although numerous clinical applications exist, the adoption of these models and methods are not common, hindered by the specialist knowledge requirements for the execution of workflows and interpretation of results. A specific focus should be increasing the access and useability for non-expert users. This should be accompanied by information regarding appropriate data collection to ensure quality simulation results. Simulation result reports, which give the clinician and patient easy to interpret answers to clinically relevant questions, should be made a priority if these modelling and simulation methods are to be widely adopted. Further integration of MSK modelling and simulations as a clinical tool cannot be achieved by MSK modelling and simulation experts alone, but requires an effort from an interdisciplinary team including modelling and simulation experts and clinicians in order to adequately balance the clinical expectations with the constraints of the modelling and simulation approaches and clinically feasible data.

**Author Contributions:** All authors (B.A.K., A.F., F.D.G. and I.J.) contributed to the conceptualization and structuring of the work, writing of the original draft, and writing and review of the revised manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Perspective* **A Conceptual Blueprint for Making Neuromusculoskeletal Models Clinically Useful**

**Benjamin J. Fregly**

**Citation:** Fregly, B.J. A Conceptual Blueprint for Making Neuromusculoskeletal Models Clinically Useful. *Appl. Sci.* **2021**, *11*, 2037. https://doi.org/10.3390/ app11052037

Academic Editor: Carlo Albino Frigo

Received: 30 December 2020 Accepted: 20 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Department of Mechanical Engineering, Rice University, Houston, TX 77005, USA; fregly@rice.edu; Tel.: +1-713-348-3212

**Abstract:** The ultimate goal of most neuromusculoskeletal modeling research is to improve the treatment of movement impairments. However, even though neuromusculoskeletal models have become more realistic anatomically, physiologically, and neurologically over the past 25 years, they have yet to make a positive impact on the design of clinical treatments for movement impairments. Such impairments are caused by common conditions such as stroke, osteoarthritis, Parkinson's disease, spinal cord injury, cerebral palsy, limb amputation, and even cancer. The lack of clinical impact is somewhat surprising given that comparable computational technology has transformed the design of airplanes, automobiles, and other commercial products over the same time period. This paper provides the author's personal perspective for how neuromusculoskeletal models can become clinically useful. First, the paper motivates the potential value of neuromusculoskeletal models for clinical treatment design. Next, it highlights five challenges to achieving clinical utility and provides suggestions for how to overcome them. After that, it describes clinical, technical, collaboration, and practical needs that must be addressed for neuromusculoskeletal models to fulfill their clinical potential, along with recommendations for meeting them. Finally, it discusses how more complex modeling and experimental methods could enhance neuromusculoskeletal model fidelity, personalization, and utilization. The author hopes that these ideas will provide a conceptual blueprint that will help the neuromusculoskeletal modeling research community work toward clinical utility.

**Keywords:** neuromusculoskeletal model; musculoskeletal model; computational model; model personalization; treatment optimization; treatment design; predictive simulation; human movement

#### **1. Introduction**

Engineers love a challenging problem, and some of the most challenging problems that engineers seek to tackle involve human health. From an engineering perspective, the human neuromusculoskeletal system has the appearance of a well-designed mechanical system. It possesses rigid segments connected by extremely low-friction rotational joints. The joints are controlled by highly efficient linear actuators whose masses are strategically located as far inboard as possible. The control system for the actuators can stabilize a multi-segment inverted pendulum flawlessly under a wide variety of challenging movement conditions. Because of these similarities to an engineered system, the human neuromusculoskeletal system lends itself well to physics-based engineering analysis and simulation. This observation raises the possibility that the same computational technologies that have revolutionized the design of airplanes and automobiles over the past 25 years could also be used to revolutionize the design of surgical and rehabilitation treatments for movement impairments. Clinical conditions that could potentially benefit include stroke, osteoarthritis, Parkinson's disease, spinal cord injury, cerebral palsy, limb amputation, and orthopedic cancer.

Though computational models representing different aspects of the human neuromusculoskeletal system have been developed for decades [1–16], such models have yet to make a positive impact on the design of treatments for movement impairments. The

neuromusculoskeletal modeling literature is full of journal articles with titles like "A model of x for doing y," but actual clinical application of the model is rarely reported. Furthermore, few, if any, details are provided on how exactly the computational model would be used to improve clinical treatment design. Nonetheless, a few research groups have generated and validated initial clinical predictions using computational models, such as how an individual with medial knee osteoarthritis should walk differently to relieve pain and slow progression of the disease [17] and whether a child with an equinus deformity due to cerebral palsy should receive gastrocnemius lengthening surgery to improve walking function [18]. Unfortunately, such examples are the exception rather than the rule. As recently stated by the Stanford OpenSim team, "The potential to use subject-specific simulations to understand the causes of movement deviations and to assess treatment options is exciting, but has not been fully realized" [19].

Realization of this untapped potential could initiate a paradigm shift in the way treatments are designed for movement impairments (Figure 1). The current treatment design paradigm has not progressed substantially beyond *generic interventions* identified by costly and time-consuming clinical trials and selected using *subjective clinical judgment*. A clinical trial collects a large quantity of movement data from a large heterogeneous patient population and uses statistical models to test whether a particular intervention induces the desired outcome in general. In contrast, a future paradigm could use *personalized interventions* identified and selected using inexpensive and efficient *objective computational models* [20,21]. This alternate approach would collect a small quantity of movement data from a specific patient and use a patient-specific computational model to test how a wide variety of possible interventions would affect the patient's functional outcome. In essence, this new approach would systematically apply virtual treatments to a virtual patient with the goal of identifying the most effective treatment for that patient. Computational treatment design using patient-specific neuromusculoskeletal models would be consistent with current clinical emphases on both personalized medicine and evidence-based medicine.

A computational approach to treatment design using personalized neuromusculoskeletal models possesses several advantages over the current treatment design paradigm. First, a computational approach would increase objectivity in treatment planning. Second, it would allow exploration of cause-and-effect relationships that map treatment decisions directly onto the patient's functional outcome and do so without burning bridges. Third, it would facilitate identification of previously unknown treatments and sensitive treatment parameters. In addition, by using physics-based rather than machine learning models, computational treatment design would extrapolate well to new situations (e.g., post-treatment patient function) for which no experimental data are available. These advantages go beyond those provided by generic neuromusculoskeletal models, which can be used to elucidate factors that are, in general, likely to contribute to the development and progression of disease and thus may serve as targets for intervention (e.g., [22,23]).

This article provides the author's personal perspective—a conceptual blueprint—for moving computational neuromusculoskeletal models to the point of clinical utility. The focus of the blueprint is on how the neuromusculoskeletal modeling research community can use commonly-employed modeling methods (i.e., rigid body dynamic, geometric, and lumped parameter) combined with commonly-available experimental data (i.e., video motion capture, ground reaction force, and electromyographic (EMG)) to achieve clinical utility. The article is structured as follows. In Section 2, motivation is provided for why computational neuromusculoskeletal models have strong potential to become clinically useful tools. Section 3 discusses challenges to progress that must be overcome if neuromusculoskeletal models are to become clinically useful and provides suggestions for how to overcome them. Section 4 provides a description of the clinical, technical, collaboration, and practical needs that must be satisfied for neuromusculoskeletal models to fulfill their clinical potential, along with recommendations for meeting them. Section 5 discusses opportunities involving more complex modeling and experimental methods that could enhance neuromusculoskeletal model fidelity, personalization, and utilization. Finally, Section 6 provides concluding thoughts on getting neuromusculoskeletal models across the threshold of clinical utility. While the emphasis of this article is on the computational design of treatments for movement impairments affecting gait (the author's primary research focus), the concepts presented should apply equally well to research involving computational design of treatments for movement impairments affecting the upper extremities and spine.

**Figure 1.** Comparison of (**a**) current and (**b**) future treatment design paradigms. The current paradigm relies on an implicit mental model in the mind of the clinician. Given clinical and imaging data and proposed treatment parameters as inputs, the implicit mental model produces a subjective prediction of post-treatment function, so different clinicians may propose extremely different treatment designs. The future paradigm replaces the implicit mental model with an explicit computational model that obeys laws of physics and principles of physiology. With this approach, movement data are added to the inputs, the explicit computational model produces an objective prediction of post-treatment function, and the entire process is wrapped in numerical optimization to identify the treatment design that will maximize the patient's functional outcome.

#### **2. Motivation for Modeling**

For neuromusculoskeletal modeling researchers, the motivation for modeling is clear there is significant potential for personalized models to improve the design of clinical treatments for movement impairments. However, clinical and basic science researchers often do not share the same perspective and are frequently skeptical about the potential benefits of computational modeling for this clinical area. This perspective is understandable given that computational treatment design for movement impairments has not yet produced numerous clinical "wins." However, this perspective also introduces funding challenges that make rapid research progress difficult. For example, in the United States, the National Science Foundation primarily uses engineering reviewers, funds innovative modeling research, does not fund clinical research, and typically funds at a level that allows limited research progress. In contrast, the National Institutes of Health uses a significant number of clinical and basic science reviewers, primarily funds clinical and basic science research, is less likely to fund innovative modeling research, but funds at a level that allows rapid research progress. Thus, a potential funding mechanism gap exists for large innovative neuromusculoskeletal modeling research projects. If more rapid advances are to be made, modeling researchers need to provide clinical and basic science reviewers with strong,

well-reasoned arguments for why model-based treatment design is worth pursuing. Below are several observations to help with building such arguments.

To begin, consider the following question: What would it take to make neuromusculoskeletal modeling a clinically useful tool for improving the design of treatments for movement impairments? To develop a computational process for treatment design, much less than one might expect. Over the past 25 years, neuromusculoskeletal modeling research has crept to the top of a precipice overlooking the "valley of death," the proverbial chasm between fundamental research on one side and clinical applicability on the other [24]. With improvements in computational modeling capabilities and computer speed, the ability to personalize a neuromusculoskeletal model to represent a specific patient and then use that model to predict an optimal treatment for that patient is no longer science fiction [17,25–27]. "We have the technology." [28] Instead, the main problem is making the computational technology readily accessible by neuromusculoskeletal modeling researchers working in collaboration with clinicians, which in turn is a problem tied to the funding mechanism gap noted above.

Is it worth trying to cross the "valley of death" in the hope that computational modeling can transform the design of treatments for movement impairments? The author argues that it is worth trying for at least four reasons. First, it is worth trying due to the significance of the clinical problems being addressed. Osteoarthritis, stroke, spinal cord injury, traumatic brain injury, and amputation affect roughly 19% of the U.S. adult population [29], with osteoarthritis and stroke being leading causes of serious long-term disability in adults worldwide [29–31]. Along with other conditions such as cerebral palsy, Parkinson's disease, and orthopedic cancer, these conditions often significantly impair movement, resulting in substantial societal costs (e.g., health care, lost productivity), an increased risk of serious secondary health conditions (e.g., heart disease, diabetes), a reduction or even loss of independence, and a decreased quality of life [29,32]. Despite the significance of these clinical problems and the uniqueness of each patient's clinical situation, treatment design has not progressed substantially beyond trial-and-error implementation of generic interventions selected through subjective clinical judgment. Even when the latest technology is used (e.g., rehabilitation robots, functional electrical stimulation), interventions are typically standardized rather than customized to the unique needs of the patient [33,34]. This state of affairs has contributed to suboptimal functional outcomes, with patients often recovering less function than desired. For example, for individuals with physician-diagnosed arthritis, approximately 44% experience significant functional limitations due to the disease [35]. Of those who receive total knee arthroplasty for knee osteoarthritis, as many as 17% are not satisfied with the outcome [36]. For individuals who suffer a stroke, only 65% regain ambulatory function, but their gait is typically slow, asymmetrical, and metabolically inefficient [37–40]. Furthermore, only 50% of individuals with walking dysfunction following stroke respond to intervention [41]. For individuals who suffer limb amputation, between 40% and 60% are not satisfied with their prosthesis [42]. For individuals who receive limb-sparing hemipelvectomy surgery for pelvic sarcoma, few are able to regain normal walking function [43,44]. We can and should do better. If individuals affected by impaired movement are to recover the most function possible, a new approach—one that is more effective, objective, and personalized—is needed for neurorehabilitation, surgical, and prosthetic treatment design.

Second, it is worth trying since the same computational design approach has revolutionized the design of products in the aerospace, automotive, heavy equipment, medical device, and numerous other industries [45–50]. In the early 1990s, Boeing was faced with the challenge of bringing a new long-range commercial airplane to market with significantly reduced development time and cost. That airplane was the Boeing 777, the most successful wide-body commercial airplane in aviation history [51,52]. To achieve its design goals, Boeing made the bold decision to reject the costly and time-consuming traditional design process involving trial-and-error design iterations performed experimentally on physical prototypes. Instead, it followed a cheaper and faster emerging design process

involving systematic design iterations performed computationally on virtual prototypes. The Boeing 777 became the first airplane designed entirely on the computer [52,53], with the final design being developed in reduced time, at reduced cost, and with increased flight range, fuel efficiency, and mechanical reliability [52,54].

Twenty-five years later, it is impossible to fly in a new commercial aircraft or drive in a new car that was not designed using computational modeling. Why? Because computational models identify better final designs than can be found through physical intuition and manual iteration [45,55,56]. As noted in a recent BBC Future news story on airplane design, "By trawling through an exhaustive set of [design] options, computers typically find ones that a human would have missed" [45]. The Food and Drug Administration (FDA) in the United States has also acknowledged the potential benefits of computational modeling for the design of medical devices and interventions [57,58]. As recently stated by Dr. Tina Morrison, Deputy Director of the FDA's Center for Devices and Radiological Health:

*"FDA recognizes the public health benefits offered by modeling and simulation, including those in the area of in silico clinical trials (using individualized computer simulation in development and or regulatory evaluation of medical products, medical devices, or medical interventions)."* [59]

Third, it is worth trying since the same computational design approach has already demonstrated success for designing personalized interventions for other clinical problems such as craniofacial bony defects [60–62] and occluded coronary arteries [57,63–65]. In both cases, the interventions are personalized using patient-specific physics-based computer models whose shapes are optimized. Unlike "black-box" machine learning models, physics-based models extrapolate and can accurately predict function for situations (such as post-treatment) outside the training boundaries. For craniofacial surgery, repeated finite element simulations allow clinician-engineer teams to design and implant personalized craniofacial prostheses that restore a patient's jaw function and aesthetics [60–62]. For coronary artery stenting procedures, large numbers of computational fluid dynamics (CFD) simulations performed by the company HeartFlow allow clinician-engineer teams to "identify significant coronary artery disease and determine the optimal treatment pathway," [63] reducing invasive coronary angiograms by 61%, unnecessary angiograms by 81%, and healthcare costs by 26%, all while improving patient quality of life and satisfaction [63,65,66]. Two related interventions designed using patient-specific physicsbased models are not far behind—nasal airway obstruction surgery planned using CFD simulations [67] and atrial fibrillation ablation surgery planned using computational electrophysiological simulations [64,68].

Fourth, it is worth trying since the same computational approach has already demonstrated success for designing at least one personalized intervention for a movement impairment—a novel rehabilitation treatment for medial knee osteoarthritis [69]. The computational design process optimized the motion of a patient-specific full-body dynamic skeletal model to predict a personalized gait modification that would minimize an external indicator of medial knee contact force. The intervention was proven to work effectively on the patient for whom it was designed, providing benefits similar to those of high tibial osteotomy surgery simply by learning to walk differently [17]. It also reduced medial knee contact force significantly in a subject implanted with an instrumented tibial prosthesis [70]. This novel intervention, which is now being investigated by research teams around the world [71–77], was identified not by subjective clinical judgment but rather by an objective computational model.

These four reasons suggest that it is both worthwhile and strategic to "help [neuromusculoskeletal modeling] researchers engage in clinical research—and cross the valley of death" [24]. However, providing skeptical grant reviewers with a strong motivation is not enough. Significant challenges to progress and significant needs must be addressed to move neuromusculoskeletal modeling to the point of clinical utility.

#### **3. Challenges to Progress**

At least five challenges help explain why neuromusculoskeletal modeling has not yet progressed to the point of clinical utility.

#### *3.1. Movement Data Alone Do Not Provide the Answer*

Back in 1932, Schwartz and Heath stated in the *Journal of Bone and Joint Surgery (American)*: "Empiricism fostered by trial and error, must continue to govern the therapy of abnormal function until measurement in some form improves the treatment of disabilities affecting the back and lower extremities" [78]. Restated more succinctly, their hypothesis was, "If we can measure it, we can fix it." How has that hypothesis worked out over the past 90 years? We can now measure more movement-related quantities than ever before. Motion can be measured routinely using video-based motion capture technology, inertial measurement units, and even raw video images. Forces exerted on the ground and the environment can be measured using various types of load cells. Pressures under feet and within prosthetic sockets can be measured using various types of pressure sensors. Muscle electrical activity can be measured using surface and fine-wire EMG sensors. Highly accurate bone motion can be measured using single-plane and bi-plane dynamic x-ray techniques. Despite these advanced measurement capabilities, the resulting movement data have yet to produce transformative improvements in the treatment of movement impairments. The one exception is surgical treatment design for gait deficits caused by cerebral palsy, where gait lab measurements often change surgical decisions [79,80]. However, even for that clinical application, different clinicians will make different surgical recommendations given the same gait data for a specific patient [81]. Thus, as noted by clinician and biomechanics researcher Richard Brand, "Although locomotion analysis has always held great promise to aid the clinician, it has never lived up to its promise" [82].

Why has the ability to measure multiple aspects of human movement not lived up to its promise? The problem is that movement data tell us what a patient does but not necessarily how to fix the patient's problem. Because the human neuromusculoskeletal system is a highly nonlinear dynamic system, it is difficult to predict intuitively how a change in system properties (e.g., changing the geometry of a bone, the length of a muscle, the attachment point of a muscle, the control strategy of the central nervous system, or the design of an assistive device) will change the patient's movement function. The relationship between muscle excitation (as measured by EMG) and muscle force is described by nonlinear activation and contraction dynamics and nonlinear muscle force-length and force-velocity relationships [83]. The relationship between muscle force and muscle moments is described by nonlinear muscle moment arms [6]. Finally, the relationship between muscle moments and resulting body motion is described by nonlinear skeletal dynamics. Each of these nonlinearities contributes to the complexity of the human neuromusculoskeletal system, making intuition about how a particular treatment will affect the patient's functional outcome unreliable. A more objective and reliable approach is needed to replace intuition.

#### *3.2. Every Patient Is Unique*

Each individual suffering from a movement impairment is unique neurologically, physiologically, and/or anatomically. Depending on the type and cause of movement impairment, individual differences can significantly impact the effectiveness of any proposed surgical or rehabilitation intervention. For example, "At present, the stroke rehabilitation field faces the challenge to tailor evidence-based treatment strategies to the needs of the individual stroke patient" [84]. For this reason, personalized rather than generic neuromusculoskeletal models are needed to support model-based intervention design for clinical conditions where patients exhibit significant heterogeneity [21,85]. Furthermore, without appropriate personalization, neuromusculoskeletal models do not reliably predict internal muscle and joint contact forces, body motion, or metabolic cost [25,86–94]. Unfortunately, the model personalization process is challenging, as some model parameter values cannot be measured directly and are only weakly observable. To make matters worse, individuals

with movement impairments are likely to have important model parameter values that differ substantially from typical values obtained from healthy individuals. As stated in more general terms by the National Academy of Engineering, "Doctors have long known that people differ in susceptibility to disease and response to medicines. But with little guidance for understanding and adjusting to individual differences, treatments developed have generally been standardized for the many, rather than the few." [95]

#### *3.3. People Change over Time*

By definition, "biomechanics" is the application of engineering mechanics to biological systems. One of the primary differences between a biological system and a man-made system is that biological systems have the capacity to adapt, either positively or negatively, over time. Such adaptation can influence both progression of and recovery from disease. However, how the human neuromusculoskeletal system will adapt over time to a specific clinical situation or intervention is difficult to simulate or predict. Apart from simulation studies of bone [96–101] and muscle [102] adaptation in response to mechanical loading, computational simulation of adaptation in the human neuromusculoskeletal system remains understudied. As noted by Richard Brand, "[D]espite powerful investigational tools, I would argue biomechanics has made a relatively minor impact in clinical practice primarily because most studies fail to account for the major distinction between living and nonliving systems: adaptability. While any study requires a clear question or hypothesis or goal, without accounting for adaptability, these studies might well be termed 'necromechanical'" [103]. Although treatments for movement impairments by definition seek to elicit a change in function over time, most biomechanical research investigates only a snapshot in time and thus might be better labeled as necromechanical.

At least two methods exist for simulating neuromusculoskeletal adaptation. The first method is harder and involves simulating the process of adaptation. This method requires the development of fundamental adaptation principles that describe how the neuromusculoskeletal system changes gradually in response to an intervention, a change in mechanical environment, or simply the progression of time. For example, Reinkensmeyer and colleagues have suggested that the ability to model time-dependent mechanisms of neuroplasticity and motor learning, which reflect neural adaptation, could help predict a patient's potential for recovery [104]. The second method is easier and involves predicting the outcome of adaptation. This method relies on identification of optimization principles that successfully predict a patient's functional outcome once recovery has plateaued. For example, Meyer and colleagues have shown that minimization of changes in the modular control signals of an individual post-stroke can successfully predict how the individual will walk under new conditions [25]. For both methods, significant additional research is needed to identify appropriate methods for modeling and simulating adaptation so that functional outcome can be predicted reliably for individual patients.

#### *3.4. Validation Is Often Weak*

Although a computational neuromusculoskeletal model is only as useful as its ability to predict reality, the research community has made only limited progress at validating its predictions of human movement and the internal forces experienced by muscles, ligaments, and articular surfaces [105]. Since muscles are the actuators of human movement, unvalidated muscle force predictions have, in turn, limited progress in using computational models to predict post-treatment patient function for different treatment scenarios. Prediction of the internal forces experienced by muscles, ligaments, and articular surfaces is challenging due to the muscle redundancy problem [106–109]. Because the human body possesses approximately three times more muscles than degrees of freedom in the skeleton, no unique mathematical solution exists for the muscle forces generated during human movement.

The main reason for limited model validation is the lack of direct in vivo measurements of internal forces before treatment and of internal forces and patient function after

treatment. In vivo measurement of hip and knee contact forces during walking and other activities has become possible through the use of instrumented implants [110–122]. Though limited to a small number of patients with implanted rather than healthy joints, these measurements provide valuable information about the in vivo loads experienced by articular contact surfaces during various functional activities. Several research groups have shared extensive human movement data sets that include in vivo hip or knee contact force measurements along with video motion capture, ground reaction force, EMG, single- or bi-plane fluoroscopy, and CT and/or MR data [121–123]. Some research groups are making extensive use of these unique shared data sets for model validation purposes. However, other research groups have not used them and consequently are now receiving negative grant proposal reviews in the United States from reviewers who are aware of these validation data sets. Direct measurement of muscle or tendon forces in vivo has historically been possible only under special circumstances [124–129]. However, researchers have recently developed a novel shear wave tensiometer to measure tendon forces in vivo during human movement [130,131]. These measurements are likely to become the gold standard for validating muscle force predictions generated by computational neuromusculoskeletal models.

An associated reason for limited model validation is that even when unique experimental data sets exist, some research groups are unwilling to share them with the research community. The main reason is the fear of generating research competition that could result in getting "scooped" on grant proposals or journal publications. To facilitate the "Grand Challenge Competition to Predict In Vivo Knee Loads," the author and colleagues [121] freely disseminated the most widely utilized human movement data sets to date with over 18,000 downloads from SimTK [132]. Since the release of the first data set in 2011, the author's research group has yet to be "scooped" in any grant proposal or journal publication, even though the data sets continue to be used extensively by the research community. While freely distributing these unique data sets has been nothing but beneficial for the author's research program, taking the same step will likely require a "leap of faith" for some research groups. Regardless, funding agencies and publishers are providing a strong nudge by requiring that collected data sets used in publications be made freely available to the entire research community [133].

To date, no shared human movement data sets have been published that provide extensive measurements of patient function both before and after treatment, along with a quantitative and qualitative description of the treatment decisions implemented for the patient. Such data sets are needed so that model predictions of a patient's post-treatment function can be validated given the patient's pre-treatment movement data and the treatment plan implemented clinically. Since inaccurate musculoskeletal model predictions could lead to the design of ineffective or even harmful treatments, it is imperative that the neuromusculoskeletal modeling research community makes model validation efforts a higher priority.

#### *3.5. Prediction of Post-Treatment Function Is Difficult*

Even given a personalized neuromusculoskeletal model, it remains challenging to predict computationally how a particular patient will function following a planned intervention. This difficulty arises for at least two reasons. The first reason involves challenges in generating predictive simulations of human movement. In 2001, Anderson and Pandy [134] published a seminal paper on predicting human walking using a full-body neuromusculoskeletal model. The walking model was three-dimensional (3D), possessed 23 degrees of freedom controlled by 54 muscles, and included deformable foot-ground contact. 10,000 h of CPU time on a cluster was required to solve the nonlinear parameter optimization problem used to predict the walking motion. That study became the standard against which all other human walking predictions would be evaluated, and it inspired subsequent efforts to generate 3D walking predictions [17,25,26,135–143]. However, even today, few research groups possess the knowledge and technical expertise needed to generate complex muscle-actuated 3D predictions of human movement.

The second reason why prediction of patient function under novel conditions is difficult involves challenges in modeling a patient's neural control strategy. To generate reliable predictions of a patient's post-treatment movement function, researchers must be able to predict how a patient's neural control strategy will change for the post-treatment conditions. When pain or fear of movement is a significant contributing factor to functional impairment, prediction of the patient's post-treatment neural control strategy is especially challenging and will require further research. Recent studies have shown that to predict patient function under novel conditions, personalized neural control models are likely to be necessary [25], and furthermore, that pre-treatment muscle synergies may facilitate the prediction of posttreatment muscle excitations [144]. While numerous studies exist that describe complex, realistic sensorimotor control models incorporating elements such as supraspinal control, modularity, central pattern generators, and proprioceptive feedback [22,145–159], these models are generic rather than personalized and are typically applied to simplified planar dynamic systems. More realistic musculoskeletal models have been coupled with simpler neural control models employing modular control [25–27,136,137,160–163] or proprioceptive feedback [164–166], but these models have yet to be evaluated in clinical treatment design scenarios. Since human movement impairments often involve a significant neurological component (e.g., stroke, spinal cord injury, cerebral palsy, even osteoarthritis), modeling approaches that can reliably predict a patient's neural control strategy following a planned intervention are important to develop.

#### **4. Description of Needs**

To address these five challenges and turn neuromusculoskeletal modeling into a clinically useful tool, the neuromusculoskeletal modeling research community must address at least four categories of needs: (1) clinical, (2) technical, (3) collaboration, and (4) practical. Clinical needs address the scope of clinical problems that can potentially be tackled with computational models and the paradigm needed for using models to design personalized interventions. Technical needs address the scope of technical developments required for model-based treatment design to become broadly feasible for the neuromusculoskeletal modeling research community. Collaboration needs address the scope of collaboration between clinicians and neuromusculoskeletal modeling researchers needed to move modeldesigned treatments from benchtop to bedside. Finally, practical needs address the scope of additional issues that must be addressed for neuromusculoskeletal modeling to become clinically useful.

#### *4.1. Clinical Needs*

Only certain clinical problems involving impaired movement are likely to benefit from a treatment design approach based on computational modeling. For example, there would be little benefit to performing computational modeling for clinical problems where "cookie-cutter" treatment plans already work well. At the other extreme, there would also be little benefit for clinical problems where genetic and molecular factors play a predominant role. Clinical problems involving orthopedic surgery, physical rehabilitation, or neurorehabilitation are likely to benefit the most from a computational approach to treatment design.

Clinical situations where a patient could benefit from computational treatment design can be identified by the presence of three elements. The first element is a "clinical measure" [82]. This element is logical since computational models produce quantitative predictions. As stated by Richard Brand, "The 'acid tests' of clinical usefulness of any measure ... are whether that measure predicts a different outcome than would be predicted without the measure, or whether the measure suggests a different treatment (or a different implementation of the same treatment) than would be recommended without the measure" [82]. In essence, an appropriate clinical measure will define a quantifiable "target" that the clinician wants to "hit" through treatment. Clinical measures that meet these criteria will be highly correlated with some direct measure of clinical outcome, such

as slowing disease progression, decreasing pain, or increasing function. Examples of clinical measures include medial knee contact force or the peak knee adduction moment for individuals with medial knee osteoarthritis or walking speed and bilateral symmetry for individuals with stroke or Parkinson's disease. Identification of the best clinical measure for any particular patient and movement impairment will require the combined input of clinicians, the patient, and researchers.

The second element is "clinical requirements." For computational models to be used for orthopedic surgery, physical rehabilitation, or neurorehabilitation treatment design, five clinical requirements should be met. First, the movement impairment should be influenced heavily by patient-specific factors. Such factors include unique neurological, physiological, or anatomical characteristics, as noted above. Second, standard treatment approaches should be either ineffective or inconsistent at restoring movement function. If standard treatment approaches already work well, then a computational approach to treatment design is unnecessary. Third, a clinical measure can be identified that can quantify initial clinical state and final treatment outcome. Without such a measure, the computational model would have no quantitative "target" to be "hit" through treatment. Fourth, clinical treatment parameters can be identified that can be changed by the clinician, can be modified in the model, and can affect the clinical measure. Identification of appropriate clinical treatment parameters requires close collaboration between clinician and researcher. Fifth, the clinical measure can be measured in the lab prior to treatment, predicted with a computational model to design the treatment, and measured again in the lab following the treatment. This process of "closing the loop" is essential for evaluating whether a computational model can predict functional outcome initially for treatments not designed with the model and subsequently for treatments designed using the model.

The third element is a "clinical model." Three model requirements should be met for a computational model to be used for clinical treatment design. First, a clinical model should be predictive, capable of predicting how the clinical measure will change in response to a proposed treatment implementation. Second, a clinical model should be minimally complex, containing only those model elements necessary to predict the clinical measure reliably for an individual patient. Third, a clinical model should be patient-specific, capable of having relevant model parameter values calibrated to the patient's pre-treatment movement and other data.

A real-life example involving computational design of a personalized rehabilitation treatment for medial knee osteoarthritis demonstrates how these three elements work together [17]. The patient studied had knee pain due to medial compartment knee osteoarthritis and was seeking an effective non-surgical treatment option. The goal became to design a modified gait pattern that reduced the patient's knee pain. The selected clinical measure was the peak knee adduction moment. This quantity is highly correlated with the rate of medial knee osteoarthritis progression [167], and a low value following high tibial osteotomy surgery is associated with the best long-term functional outcome [168]. For clinical requirements, knee loading is heavily influenced by patient-specific leg alignment and walking pattern, standard non-surgical treatment approaches are ineffective at relieving pain and restoring walking function, the peak knee adduction moment can quantify initial clinical state and final treatment outcome, the subject's walking motion can be modified and will change the peak knee adduction moment, and finally the peak knee adduction moment can be measured in the lab prior to treatment, predicted with a computational walking model to design treatment, and measured again in the lab following treatment. The clinical model predicted the patient's peak knee adduction moment for different modified walking motions, was minimally complex by modeling only the skeleton (i.e., no muscle or neural control models) and omitting foot-ground contact models, and utilized a patient-specific kinematic structure and mass distribution. An optimization problem was formulated to predict the patient-specific gait modifications that would minimize the patient's knee adduction moment while still producing a normal-looking walking motion. The optimization predicted that a modified gait pattern involving knee medialization

during stance phase could produce a 34% reduction in the patient's peak knee adduction moment. After learning to perform the predicted gait modification, the patient was able to achieve a 37% reduction when retested in a gait lab [17]. This reduction was comparable to that expected from high tibial osteotomy surgery [169] and decreased the patient's risk of osteoarthritis progression by a factor of 10 [167]. The patient incorporated the modified gait pattern into his daily life, experienced significant pain reduction, and was able to return to running and recreational sports activities.

#### *4.2. Technical Needs*

While the potential is great for using computational models to design more effective treatments for movement impairments, realization of this potential is currently limited by two technical challenges. The first technical challenge involves model personalization. As noted above, for clinical conditions where every patient is unique, personalized models are needed to support model-based intervention design. For movement impairments in particular, it can be beneficial to personalize four aspects of a patient's neuromusculoskeletal model using the patient's pre-treatment movement and neurophysiological data [25,26], where the model elements requiring personalization will depend on the clinical problem at hand. These four aspects are: (1) joint model personalization, where parameters defining patient-specific joint centers and functional axes are calibrated [170–176], (2) muscle-tendon model personalization, where parameters defining patient-specific Hill-type muscle force generation [83,88,90,177–184] and (if desired) surrogate musculoskeletal geometry [185,186] are calibrated, (3) foot-ground contact model personalization (for walking impairments), where parameters defining patient-specific deformable ground contact force generation are calibrated [187], and (4) neural control model personalization, where parameters defining a patient-specific neural control structure are calibrated [25].

While computational methods already exist for personalizing each of these four aspects of a neuromusculoskeletal model (Table 1, 2nd through 5th columns), few research groups are personalizing their models. Furthermore, those that do typically personalize only one or two aspects of their models (review Table 1). The main reason is simple—model personalization is hard, requiring development of specialized optimization approaches that identify model parameter values that fit a patient's movement and neurophysiological data as closely as possible. Even within the author's research group where all four aspects are personalized [25,26], it takes a new Ph.D. student six months or longer to learn how to use the lab's existing model personalization processes. While musculoskeletal modeling software such as OpenSim [19,188] and AnyBody [189] facilitate the process of constructing models and performing standard analyses (e.g., model scaling, inverse kinematics, inverse dynamics, forward dynamics, static optimization), these tools do not currently provide standard methods for personalizing these four model aspects. Realizing this gap, the Stanford OpenSim team recently noted, "More development is needed to streamline the process of creating and validating simulations of individuals with impairments" [19].

The second technical challenge involves treatment optimization. Treatment optimization adjusts neural control signals along with surgical, neural control, internal implant, and/or external device parameters in the patient's personalized model to achieve a specified treatment goal (e.g., maximize walking speed and symmetry). As recently stated by the OpenSim team, "The prediction of outcome due to treatment or intervention (surgery, physical training, biofeedback, etc.) remains the ultimate goal of musculoskeletal modeling and simulation" [19]. However, even if an appropriately personalized neuromusculoskeletal model is available, using that model to predict how the patient will function following a specific type of intervention or a specific implementation of the intervention remains challenging. Since Anderson and Pandy's foundational predictive walking simulation [134], a number of research groups have developed the capability to predict human movement under novel conditions (Table 1, last column). Despite these advances, most predictive simulations of human movement use scaled generic two-dimensional models rather than personalized 3D models and are rarely performed for treatment optimization purposes. Furthermore, most predictive movement simulations solve a tracking optimization problem, where the goal is to reproduce an experimentally measured motion rather than predict a new motion for a post-treatment situation. Similar to model personalization, the main reason why treatment optimization has not gained widespread use is that it is hard, again requiring specialized optimization approaches, in this case to predict how the patient will function after a planned intervention is applied to the patient's model.

**Table 1.** Overview of neuromusculoskeletal modeling research involving model personalization and motion prediction using optimization methods.


Realizing this challenge, neuromusculoskeletal modeling researchers have explored different nonlinear optimization methods derived from the field of optimal control for performing predictive simulations of human movement. Early efforts employed a direct shooting method [8,246]. For this method, the optimization design variables are model controls (e.g., muscle excitations), repeated forward dynamic simulations are performed explicitly through numerical integration, and individual time frames are solved sequentially by time marching. The disadvantages of this approach are that repeated forward dynamic simulations are susceptible to numerical integration drift and other integration problems, especially when intermittent foot-ground contact is involved, plus unstable movements such as walking cannot be predicted without some form of stabilizing feedback control. To circumvent these problems, researchers have recently converged on direct collocation as the preferred method for predicting human movement [25,26,190,232,242,243,245,247,248]. For this method, the optimization design variables are model controls as well as states, repeated forward dynamic simulations are performed implicitly as part of the optimization problem formulation, and all time frames are solved simultaneously, thereby eliminating time marching. While this approach produces a larger nonlinear optimization problem, it eliminates numerical integration drift and stability problems, and it is generally more reliable for predicting new motions. While some researchers have implemented their own direct collocation optimal control methods [190,191,242,243,245,247,249–251], implementation of these methods is non-trivial, leading other researchers to use academic direct collocation software such as GPOPS-II [25,26], CasADi [135,136], and OpenSim Moco [248]. Despite the availability of these programs, incorporating a generic or personalized neuromusculoskeletal model into one of these packages remains challenging, limiting the use of this approach by the neuromusculoskeletal modeling research community.

The common thread in these two technical challenges is ease of use, which includes ease of implementation. Ideally in a matter of hours, researchers would be able to personalize a neuromusculoskeletal model to a patient's movement and neurophysiological data

and then use the model to predict the optimal treatment for the patient. Given the computational technologies and methods that exist today, the main reason this ideal situation is not already a reality (at least for some clinical problems) is the lack of the software tools necessary to make the model personalization and treatment optimization processes fast and easy to set up and perform. As stated in a recently review article on the use of musculoskeletal models in clinical practice, "A specific focus should be increasing the access and usability for non-expert users" [252]. For model personalization and treatment optimization to become widespread, the barriers to entry need to be lowered substantially and the entire process needs to be made easily accessible to the broader neuromusculoskeletal modeling research community.

#### *4.3. Collaboration Needs*

Even if the technical needs described above were satisfied, a significant collaboration challenge would still exist. This challenge involves finding clinical and neuromusculoskeletal modeling researchers who have the mindsets needed to work effectively together.

An illustration is helpful for understanding the joint mindset needed for effective collaboration. The illustration involves two children playing in the same sandbox and highlights the need for a concept called "shared intellectual investment." One child represents a clinician, the other child represents a modeler, and the sandbox represents the clinical problem at hand. In the least productive scenario, the clinician says to the modeler, "I don't need modeling capabilities. Why are you in *my* sandbox?" Similarly, the modeler says to the clinician, "I don't need clinical perspective. Why are you in *my* sandbox?" While rarely stated explicitly, this scenario is probably the most common, contributing greatly to the lack of progress. In a slightly more productive scenario, the clinician says to the modeler, "Come work on my clinical problem. Come build in *my* sandbox." In this case, intellectual investment in solving the clinical problem comes primarily from the clinical side. Similarly, the modeler says to the clinician, "Come use my models. Come build in *my* sandbox." Here, intellectual investment comes primarily from the modeling side. While this scenario is better than the first one, it is still far from ideal, as each side sees itself as independent from the other. In the most productive scenario, the clinician and modeler say to each other, "How could modeling help clinical problems? Let's build in the sandbox *together*." For this scenario, the sandbox becomes a shared space of shared intellectual investment, where the clinician wants to work closely with the modeler and vice versa. Both sides enter into the sandbox together with humility and appreciation for the capabilities and perspective that the other side brings to the problem. Though such collaborations exist today, they are far less common than they need to be for computational models to make a significant positive impact on treatment design for movement impairments.

Creation of a collaboration marked by "shared intellectual investment" requires a commitment from both the clinician and the modeler. The clinician needs to be willing to teach the modeler about the details of the clinical problem and to document and measure the treatment decisions implemented in the patient's treatment plan. That information is critical to allow the modeler to evaluate post-treatment predictions of patient function given a model of the patient constructed from pre-treatment data and the treatment decisions implemented by the clinician. On the other side, the modeler needs to be willing to teach the clinician about the quantities that models can reliably predict and to learn as much as possible about the issues present in the clinical problem. While the urgency of daily required tasks (e.g., surgery, clinical care, teaching, research, administration) makes this additional effort inconvenient for both sides, it is a long-term investment with the potential to pay high dividends.

#### *4.4. Practical Needs*

Beyond the previous three categories of needs, a fourth category involves practical considerations that are not rocket science but rather strategic decisions. One practical consideration is the focus of neuromusculoskeletal modeling research. Engineering re-

searchers love technology development. It can be difficult for us to know when to stop developing technology and start applying the technology we have been developing to important problems. It is the author's belief that a subset of movement impairment clinical problems exists today that could benefit greatly from existing computational modeling methods. A greater focus on identifying those problems, rather than on developing more technology, could lead to the initial "clinical wins" that would help propel model-based treatment design forward.

Another practical consideration is the need for extensive experimental data sets that would allow neuromusculoskeletal modeling researchers to develop and evaluate model-based treatment design methods. Such data sets would be collected from specific patients and include pre- and post-treatment (following plateau in recovery) movement and neurophysiological data, along with documented and measured treatment decisions implemented by the clinician. The specific experimental data needed would depend on the clinical problem but could include video motion capture, EMG, and (for walking) ground reaction data, along with medical imaging data (e.g., CT and/or MR). Pre- and posttreatment imaging data could be valuable for improving personalization of the patient's pre-treatment model and measuring certain treatment decisions made by the clinician (e.g., locations of bone cuts). Having pre- and post-treatment data and treatment decisions available for specific patients would allow researchers to apply the actual treatment decisions to the patient's pre-treatment model and predict the patient's post-treatment function. Differences between measured and predicted post-treatment function would then be used to improve and refine the modeling methods without the risk of negatively impacting the patient's clinical care. Once retrospective prediction of post-treatment patient function can be performed reliably, prospective prediction of post-treatment function could begin to be utilized clinically as long as any necessary regulatory approvals had been addressed. The author is aware of only one study to date that has "closed the loop" by validating a prospective patient-specific model-based prediction of post-treatment function using post-treatment movement data collected from the patient [17].

Since few neuromusculoskeletal modeling research groups possess the financial, experimental, and clinical resources needed to generate pre- and post-treatment movement and neurophysiological data, the research community will need well-curated data sets to begin learning how to "close the loop" for model-based treatment design. One way to accelerate the generation and distribution of such data sets would be to organize international competitions similar to the "Knee Grand Challenge Competition" organized previously by the author and colleagues [121]. For example, modelers and clinicians together could develop a "Cerebral Palsy Grand Challenge Competition," "Stroke Grand Challenge Competition," or "Pelvic Sarcoma Grand Challenge Competition." Competition organizers would provide pre-treatment movement and neurophysiological data collected from a specific patient, along with details of the treatment plan implemented by the clinician. Competing research teams would create a personalized neuromusculoskeletal model of the patient using the pre-treatment data, implement the patient's actual treatment plan in the patient's model, and then generate a blinded post-treatment prediction of the specified clinical measure. The blinded prediction would be submitted to the competition organizers, and only then would competitors be given the patient's post-treatment movement and neurophysiological data, including the specified clinical measure, to evaluate their predictions and identify the weak links in their modeling methods. The best blinded predictions, and the subsequent modeling modifications that improved the predictions, would be presented in a special session of an annual conference and published in a special issue of a relevant journal, similar to the approach used successfully for the "Knee Grand Challenge Competition." Following several competitions, designing and evaluating treatments for movement impairments would become common practice for at least a segment of the neuromusculoskeletal modeling research community.

#### **5. Opportunities for Enhancement**

Thus far, this article has focused on how to get neuromusculoskeletal modeling to the point of clinical utility using commonly-employed modeling methods (i.e., rigid body dynamic, geometric, and lumped parameter) and commonly-available experimental data (i.e., video motion capture, ground reaction force, and EMG). This focus is consistent with the well-known principle of minimum model complexity–engineers should use the minimum complexity model that can reliably predict the quantities of interest. In the author's opinion, significant clinical utility can be achieved in the near future without moving beyond commonly-used and minimally-complex modeling and experimental approaches. Furthermore, personalized modeling and predictive simulation capabilities that are based on commonly-used approaches are the most likely to gain widespread use by the neuromusculoskeletal modeling research community. At the same time, if researchers are willing to go beyond common approaches, then related research areas employing more complex modeling and experimental methods could enhance model fidelity, personalization, and utilization.

#### *5.1. Enhanced Model Fidelity*

Opportunities for enhancement exist for generating higher fidelity neuromusculoskeletal models that capture short- and long-term tissue-level behavior. For some clinical problems where functional outcome depends on tissue-level stresses and strains, rigid body skeletal models controlled by lumped-parameter Hill-type muscle-tendon models will not be sufficient. To capture tissue-level effects in muscles, ligaments, bones, and/or articular cartilage, finite element (FE) models with deformable tissue properties are the logical choice [102,253–263]. In muscles, for example, non-uniform tissue-level behavior could be important for predicting muscle hypertrophy or injury in response to exercise, the effects of aging and disuse, the progression of muscular dystrophy, and muscle function in microgravity. In ligaments, it could be important for predicting failure mechanisms and the best way to perform reconstructions. In bones, it could be important for predicting the risk of fracture due to unfavorable loading conditions, osteoporosis, or joint replacement. In articular cartilage, it could be important for predicting conditions that initiate, accelerate, or slow the development of osteoarthritis.

For more complex clinical problems where functional outcome depends on how muscle, bone, and/or cartilage tissue changes over time, FE models can be combined with either empirical or agent-based adaptation models. As noted earlier, bone adaptation to an altered loading environment is the classic historical example of adaptive simulations. However, recent work has extended adaptive FE simulation methods to muscle and articular cartilage as well. For example, 3D FE models of muscle have been combined with agentbased adaptation models to simulate disease-related changes in muscle function caused by muscular dystrophy [255] and to propose improved tissue regeneration approaches to treat severe muscle injuries [102]. For articular cartilage, a specimen-specific 3D FE model of a patellofemoral joint has been coupled with an empirical adaptation model to predict articular cartilage wear over time during in vitro testing [264]. The challenge with adaptive simulations is that a new FE analysis must be performed each time the adaptation model updates the FE model's tissue properties. Furthermore, since it is unrealistic to simulate thousands or millions of motion cycles, the number of cycles for which the most recent FE simulation results can be extrapolated must be determined (see [265,266] for further details).

#### *5.2. Enhanced Model Personalization*

Opportunities for enhancement also exist for personalizing important properties of musculoskeletal models. Personalization of musculoskeletal geometry can be enhanced using imaging data and recently published software tools. These tools, which include NMSBuilder [267], the MAP Client [268], and STAPLE [269], can turn patient imaging (e.g., CT, MR) and/or surface marker data into personalized geometric OpenSim models (for a comparison of the capabilities of these three tools, see [270]). Together, such tools allow for generation of personalized bone geometries, muscle lines of action, and/or joint definitions, significantly improving the speed and accuracy with which personalized geometric musculoskeletal models can be created. A recent study has also provided a detailed step-by-step procedure for creating subject-specific lower limb musculoskeletal geometric models from MR imaging data [271].

Personalization of model parameter values can be enhanced using probabilistic modeling methods. Probabilistic modeling samples combinations of model inputs to characterize the distribution of corresponding model outputs. This approach can help identify model parameter values to which the outputs of interest are sensitive, thereby revealing which parameter values are worth personalizing and which are not [272,273]. For example, one probabilistic modeling study showed that when calculating lower extremity inverse dynamic joint moments during walking, errors in body segment mass properties have little influence, while errors in joint positions and orientations can have a significant influence [91]. Other probabilistic modeling studies have investigated which parameters in Hill-type muscle-tendon models and geometric musculoskeletal models are the most important to personalize when estimating lower extremity muscle forces during walking [274–277]. These studies have reported that optimal muscle fiber length, tendon slack length, and peak isometric force in Hill-type muscle-tendon models [274–276] and the attachment points of the iliacus and psoas muscles, but not those of other muscles [277], are the most critical to calibrate to subject data when estimating muscle forces.

Personalization of model parameter values can also be enhanced using various imaging methods. In some cases, imaging methods permit direct measurement of model parameter values. For example, magnetic resonance imaging (MRI) can be used to measure in vivo muscle volumes [278], which can then be converted to muscle peak isometric forces. Diffusion-tensor MRI can be used to perform in vivo measurements of muscle pennation angle and fiber type distributions [279,280]. Ultrasonography can be used to identify subject-specific force-strain parameters for tendons and optimal fiber length and pennation angle for muscles [281–283]. In other cases, imaging methods provide novel in vivo measurements that could improve estimation of model parameter values using optimization methods. For example, tendon forces measured in vivo by shear wave tensiometry [130,131] could be employed as additional constraints when personalizing muscle-tendon model parameter values using EMG, kinematic, and kinetic data.

#### *5.3. Enhanced Model Utilization*

Finally, opportunities for enhancement exist for making utilization of complex models faster and easier through surrogate modeling methods. Surrogate modeling is a special case of supervised machine learning that fits a computationally "fast" model to the inputoutput characteristics of a computationally "slow" model, simulation, or optimization (for a broader perspective on the use of machine learning methods for personalizing neuromusculoskeletal models, see [284]). The surrogate model essentially becomes a fast black-box model that replaces the slow original model, simulation, or optimization.

The surrogate modeling process generally involves a sequence of six steps [285]. First, identify one or more outputs of interest from the model, simulation, or optimization. Second, identify the inputs that affect the outputs of interest, where the number of inputs will ideally be small (e.g., <10) so that the surrogate model fitting process will not become overly complex. Third, sample combinations of inputs (called "sample points") within pre-defined bounds using design of experiments (e.g., Latin hypercube) to spread out the sample points evenly throughout the (possibly high-dimensional) design space. Fourth, perform an analysis, simulation, or optimization for each sample point to generate the corresponding outputs of interest. Fifth, fit each output of interest as a function of the sample point inputs using a black box fitting method (e.g., polynomial response surfaces, Kriging, support vector machines, artificial neural networks). Sixth, evaluate the accuracy of the resulting surrogate models using sample points not included in the fitting process. While it

can take hours or days of computation time to evaluate all of the input combinations, this computational cost is paid only once upfront, and the final surrogate models will be orders of magnitude faster.

For neuromusculoskeletal modeling applications, surrogate models have two primary benefits. First, surrogate models can make "slow" computational models "fast" for use within other "fast" models. In the context of neuromusculoskeletal modeling, the "fast" model is typically a rigid body dynamic skeletal model controlled by lumped-parameter muscle-tendon models, while the "slow" model is typically a geometric musculoskeletal model or a tissue-level elastic continuum model. Since the mechanics of the slow and fast models interact, the two models should ideally be simulated simultaneously rather than sequentially. Surrogate models have been developed to calculate muscle-tendon lengths and moment arms produced by geometric musculoskeletal models [185,186], stresses in foot tissues produced by FE models [192], strain fields in long bones produced by FE models [286], and contact forces, stresses, and/or wear in natural and artificial knees produced by FE or elastic foundation models [220,287–292].

Second, surrogate models can make "slow" computational simulations or optimizations "fast" for clinical implementation. Performing treatment design iterations with a physics-based model can require extensive computation time. If the clinician is satisfied with the final predicted treatment design, then long computation time may not be a problem. However, if the clinician would like to interact with the model and observe in real time the effect of different treatment decisions, then long computation time will be unacceptable. To circumvent this problem, one can fit a surrogate model to the outputs (e.g., clinical measures) of repeated simulations or optimizations [173,293]. The resulting surrogate models would allow the clinician to evaluate different treatment scenarios in real time. This approach has been used successfully with FDA approval by HeartFlow [63], where the barrier to clinical implementation was excessive computation time for thousands of CFD simulations. The same surrogate models can also be used to perform ultra-fast optimizations [285]. Clinicians and engineers working together could rapidly investigate tradeoffs in functional outcomes for different decisions about the relative importance of different simulation outputs included in the cost function.

#### **6. Conclusions**

In this article, the author has provided a conceptual blueprint for how neuromusculoskeletal models can become clinically useful. At the core of this blueprint is the need to make model personalization and treatment optimization easy to implement and easy to use by the broad neuromusculoskeletal modeling research community. In addition to motivating the use of neuromusculoskeletal models to improve treatment design for movement impairments, the author has highlighted challenges to progress along with suggestions for overcoming them, a description of needs with recommendations for addressing them, and opportunities for enhancement involving more complex modeling and experimental methods. In the short term, if a subset of movement impairment problems can be identified to which existing neuromusculoskeletal modeling methods could be applied effectively, the field could generate some initial clinical "wins" that would help propel it across the threshold of clinical utility and open up a new paradigm for treatment design. In the long term, model-based treatment design for movement impairments will need the endorsement of regulatory (e.g., the FDA in the United States, the European Medicines Agency in Europe) and health technology assessment (e.g., the Health Technology Assessment Programme in the United Kingdom) agencies, which would certify the clinical utility of the approach and facilitate crossing the "valley of death." If computational modeling and simulation can do for the design of movement impairment treatments even a fraction of what they have done for the design of airplanes and automobiles, their clinical impact will be transformative.

**Funding:** The author would like to thank the Cancer Prevention and Research Institute of Texas (grant RR170026) for supporting this work.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The author thanks his clinical collaborators Carolynn Patten, Valerae Lewis, Gerard Francisco, and James Chang for valuable discussions that were instrumental in forming many of the ideas in this article. He would also like to thank Richard Brand for inspiring the author's search for a formal model-based treatment design paradigm.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **The Effects of External Loads and Muscle Forces on the Knee Joint Ligaments during Walking: A Musculoskeletal Model Study**

**Carlo Albino Frigo 1,2,\* and Lucia Donno <sup>1</sup>**


**Abstract:** A musculoskeletal model was developed to analyze the tensions of the knee joint ligaments during walking and to understand how they change with changes in the muscle forces. The model included the femur, tibia, patella and all components of cruciate and collateral ligaments, quadriceps, hamstrings and gastrocnemius muscles. Inputs to the model were the muscle forces, estimated by a static optimization approach, the external loads (ground reaction forces and moments) and the knee flexion/extension movement corresponding to natural walking. The remaining rotational and translational movements were obtained as a result of the dynamic equilibrium of forces. The validation of the model was done by comparing our results with literature data. Several simulations were carried out by sequentially removing the forces of the different muscle groups. Deactivation of the quadriceps produced a decrease of tension in the anterior cruciate ligament (ACL) and an increase in the posterior cruciate ligament (PCL). By removing the hamstrings, the tension of ACL increased at the late swing phase, while the PCL force dropped to zero. Specific effects were observed also at the medial and lateral collateral ligaments. The removal of gastrocnemius muscles produced an increase of tension only on PCL and lateral collateral ligaments. These results demonstrate how musculoskeletal models can contribute to knowledge about complex biomechanical systems as the knee joint.

**Keywords:** knee joint ligaments; walking; dynamic simulation; musculoskeletal modelling; knee joint biomechanics

#### **1. Introduction**

The knee joint is one of the most complex articulations in our body. During movement, the relative displacements of distal femur and proximal tibia are determined by the shape of bone surfaces and by ligaments and muscle forces. The femoral condyles and the tibial plateau are not at all congruent and allow the proximal tibia to slide forward and backwards and rotate around its longitudinal axis (internal/external rotation). These movements usually occur during flexion/extension. Other movements like mediolateral displacement and adduction/abduction are usually prevented by proper tension of ligaments. During walking the whole structure is subjected to external forces due to foot-ground interaction, gravity and inertia forces [1,2]. The internal forces required to counteract the external ones are relatively high. For example, the tibiofemoral contact force is in the order of 200–300% of body weight during walking (forces measured through instrumented joint prostheses [3]) where the muscle and ligament forces contribute to them for approximately 1/3 and 2/3, respectively (contributions estimated through models [4]). Despite the demand for a high load-bearing capacity, the knee joint allows a considerable range of motion. That means that the ligaments must be stiff enough to prevent joint dislocation under load, but they are positioned in such a way that the main flexion-extension movement is not restricted.

**Citation:** Frigo, C.A.; Donno, L. The Effects of External Loads and Muscle Forces on the Knee Joint Ligaments during Walking: A Musculoskeletal Model Study. *Appl. Sci.* **2021**, *11*, 2356. https://doi.org/10.3390/ app11052356

Academic Editors: Wanda Lattanzi and Marco Parente

Received: 26 January 2021 Accepted: 2 March 2021 Published: 6 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The analysis of such a complex structure has been carried out by many authors through different approaches. Studies on cadaver specimens [5–7] have the advantage of dealing with the real anatomy and the real mechanical properties of tissues. However, they are difficult to realize and require special equipment to reproduce a physiological movement and to apply loads and muscle forces. On the other hand, another important approach, the "finite elements" method [8–10], could provide the estimation of all the structural variables in great detail, but the many uncertainties about the biological material properties and the need for an arbitrary definition of the boundary conditions pose a problem of results reliability [9,11]. Direct quantification of in vivo knee motion can be achieved through special techniques that extract the pose of the bones from single-plane fluoroscopy [12,13]. This technique is particularly effective for analyzing the motion of implanted knee joint prostheses [14], but its application is limited to a narrow volume of acquisition. In vivo measurement of inter-joint contact forces during locomotion and other motor acts have been obtained by implanting special instrumented joint prostheses [15,16]. The obtained results are considered as the reference for subsequent study validations [17]. Other interesting quantities like the tension of the ligaments cannot be obtained by this method. Attempts to measure them in vivo [18,19] have been made, but many technical and ethical problems occur with this method, and practically the only ligament deformations can be obtained. Musculoskeletal models have the potential to permit the estimation of forces internal to the joint, provided that a suitable representation of the joint structure is performed. Pioneering studies have been made by Morrison et al. [20] who analyzed the forces in the knee during walking. The knee was modelled as a single axis joint and in each instant along the walking cycle only a reduced number of muscles, ligaments and contact forces was assumed to be active, so that the number of unknowns was equal to the number of equilibrium equations. Collins and O'Connor [21] considered the possibility of antagonistic muscles co-contractions while computing the internal knee joint forces during walking, but their model was restrained to the sagittal plane. More recently, several publications have faced the problem with a comprehensive modelling approach [4,6,22–28]. In general, musculoskeletal models treat the bones as rigid bodies and represent the viscoelastic tissues as a combination of spring-like elements [10,28–32] with non-linear characteristics. These are dramatic simplifications, but from another point of view, they offer the possibility to simulate complex dynamic situations without excessive requirements of computational resources and time. The implementation of a musculoskeletal model of the knee joint is not so straightforward, since many details are to be defined that can strongly affect the results. First, the anatomical structure has to be designed, and this can be done based on medical images (MRI); second, the mechanical properties of the soft tissues, ligaments and cartilage have to be quantified; then, the kinematic conditions, the external loads (ground reactions and inertia forces and moments) and the muscle forces applied to the system have to be determined. Finally, the model has to be animated by the input data so that the internal loads (ligaments and contact forces) are computed based on the dynamic equilibrium principles. In the present work, the procedure has been implemented starting from a previously developed model [29], which has been improved in terms of number and position of the ligaments, adaptation to natural gait and determination of muscle forces. The objective of this work was to analyze the ligament tensions during walking and to understand how they are affected by muscle contractions. The problem is of clinical relevance for many situations in which the knee ligaments are injured, or they must be sacrificed because of knee joint arthroplasty. Our findings could help identify compensatory attitudes to be adopted to prevent ligament overloads or identify which muscles should be used to compensate for surgical removal of a ligament. Knowing the effect of muscle contractions on ligament tensions may suggest proper rehabilitation programs aimed at training appropriate muscle control.

#### **2. Materials and Methods**

#### *2.1. The Model Implementation*

Our model was implemented on the SimWise-4D platform (Design Simulation Technologies, DST, Canton, MI, USA). This software is designed for the analysis of complex dynamic systems, in which rigid bodies are connected by different kinds of joint constraints, and the interaction between object surfaces is solved in the modality called "collision". By this software, we have created a knee joint model which includes bones, ligaments and muscles. As already described in a previous paper [29], models of femur, tibia and patella were obtained through segmentation of MRI using commercial software (Amira 5.3.3, Visage Imaging, Inc., San Diego, CA, USA). MRI images referred to a Caucasian male (age: 42 years, body height: 1.72 m, body mass: 70 kg), with no signs of pathology. Contact surfaces were smoothed by using Meshmixer (version 3.5.47, Autodesk, San Rafael, CA, USA). To reduce the simulation time, portions of the distal femur and proximal tibia were split from the rest of the bones and then attached to them by rigid constraints (Figure 1). In this way, the algorithms for surfaces interaction had to work on these small bone portions only. The friction coefficient between tibia and femoral surfaces was set at 0.01 [33]. Ligaments were represented by non-linear springs attached to the bones at specific points that were identified with the support of a text of joint physiology [34] and other publications [6,22,26,29,35]. Several springs were used for each ligament, to represent the different fascicles whichmay be tensioned in different working conditions. The whole ligaments structure included (see Figure 1): (a) the anterior cruciate ligament (ACL) subdivided into the two compartments, anterior (ACL-a) and posterior (ACL-p); (b) the posterior cruciate ligament (PCL) subdivided in anterior (PCL-a) and posterior (PCL-p) bundles; (c) the lateral collateral ligament (LCL) subdivided into three components, LCL-1, LCL-2 and LCL-3 (respectively anterior, intermediate and posterior); (d) the medial collateral ligament (MCL), composed of the superficial bundles MCL-a, MCL-i and MCL-p (anterior, intermediate and posterior) and by the deep components MCL-da (anterior) and MCL-dp (posterior); (e) the anterolateral fibrous capsule (Cap-Ant-L); and (f) the two portions of the posterior capsule: Cap-Post-L and Cap-Post-M (lateral and medial, respectively).

**Figure 1.** The two bone portions considered for the contact algorithm: the distal femur (red, that was rigidly attached to the femur) and the proximal tibia (blue, rigidly attached to the tibia). The forces produced by the interaction of the surfaces are transmitted to the constraints that can measure the resultant force and the moment. The ligaments connecting femur and tibia are represented by straight springs having non-linear characteristics (see text).

The elastic characteristics of these springs were expressed in the form adopted by many authors [6,28,36] that included a quadratic rise of the force up to a predefined strain limit, and a linear behavior for larger deformations:

$$f = \begin{cases} \ 0.25 \text{K} \varepsilon^2 / \varepsilon\_l, & 0 \le \varepsilon \le 2\varepsilon\_l \\\ \ \text{K}(\varepsilon - \varepsilon\_l), & \varepsilon > 2\varepsilon\_l \\\ 0, & \varepsilon < 0 \end{cases} \tag{1}$$

*K* is the stiffness (N/m) in the linear section of the curve; *ε* is the strain, (L − L0)/L0; *εl* is the strain limit that corresponds to the intercept of the linear tract with the abscissa. According to [6] the strain limit was assumed *ε<sup>l</sup>* = 0.03. The separation of the quadratic behavior from the linear behavior occurs at *ε* = 2*ε<sup>l</sup>* . One fundamental parameter that strongly conditions the ligament force is the rest length L0. According to [26], we assumed that in a reference position corresponding to knee full extension the strain of the ligament *εre f* was known. These values are reported in the cited publication. As a consequence, knowing the reference length in our model (*Lre f*), the rest length *L*<sup>0</sup> was computed as

$$L\_0 = \frac{L\_{ref}}{\varepsilon\_{ref} + 1} \tag{2}$$

In our simulations, we adopted the values of *εre f* reported in [26]. The only exception was the PCL-a for which the reported value was *εre f* = −0.21, which appeared to be a mistake since the values for most of the other ligaments were in the range 0.02–0.06, in absolute. We assumed for this ligament the value *εre f* = −0.03 by assuming that the ligament was slack in the reference position and would be recruited when stretched above the 3% of the reference length. According to the mentioned authors, the MCL-da was also considered slack in the reference position, with a reference strain *εre f* = −0.08. The anterolateral capsule, not taken into account in the mentioned paper, was assumed to have a small deformation in the reference position (*εre f* = −0.002), considering that it will lengthen during knee flexion. All the strain values, reference lengths and rest lengths are reported in Table 1. The difference between the reference length *Lre f* and the rest length *L*<sup>0</sup> is reported in the fifth column. The values range from 2.67 mm (MCL-i) to −3.24 mm (MCL-da). A negative value means that the ligament was slack in the reference position.

**Table 1.** The ligaments characteristics adopted in the model: *εre f* and *Lre f* are respectively the strain derived from [26] and the ligament length in our reference position (straight knee); *L*<sup>0</sup> is the rest length of the ligament; *K* is the stiffness (= *F*/*ε* ) in the linear portion of the ligament characteristics; ∆*L* and *Fre f* are, respectively, the elongation and the ligament tension in the reference position. ACL: Anterior Cruciate Ligament, PCL: Posterior Cruciate Ligament (a: anterior bundle, p: posterior bundle); LCL: Lateral Collateral Ligament (1,2,3: anterior, middle, posterior component respectively); MCL: Medial Collateral Ligament (a: anterior, i: intermediate, p: posterior, da: deep anterior; dp: deep posterior components); CAP: Capsular Ligament (Ant-L: anterolateral; Post-L: posterolateral; Post-M: posteromedial components).


As to the *K* parameter (stiffness in the linear section of the force-strain characteristics), the values reported by [26] were adopted in our model. For those ligaments that appeared as single units in the mentioned paper but in our model were split, the stiffness of each bundle was obtained by dividing the reported value by the number of bundles constituting

the ligament. This occurred for LCL and posterior capsule. The stiffness values are reported in the sixth column of Table 1. The resulting pre-tensions of the ligaments in the reference position are reported in the last column. These forces ranged from 0 (slack ligaments: PCL-a, MCL-da) and 45 N (posterior capsule: CAP-Post-L and CAP-Post-M).

Once the ligament properties were defined, the muscles were accurately positioned in the model.

Each muscle was represented by a force actuator that could produce a controlled concentric force. As to the extensor mechanism, the rotula was reproduced by a cylinder sliding on the femoral trochlea (see Figure 2). This simplification does not affect the functioning of the knee extensor mechanism and helps reducing the computational time. The rotula was attached to the tibial tuberosity by a rope representing the patellar tendon. At the other extremity, the rotula was attached to a series of smaller cylinders representing the quadriceps tendon. These cylinders were connected to each other by spherical joints and could slide on the trochlear surface when the knee flexion was sufficiently high.

**Figure 2.** The muscles represented in the model: RF (Rectus Femoris), VI (Vastus Intermedius), VM (Vastus Medialis), VL (Vastus Lateralis). These muscles constitute the quadriceps (Quad). ST (Semitendinosus), SM (Semimembranosus), BF-lh (Biceps Femoris–long head), BF-sh (Biceps Femoris–short head: they constitute the Hamstrings (Ham). GaL (Gastrocnemius Lateralis), GaM (Gastrocnemius Medialis): they constitute the Gastrocnemius (Gas).

> Four actuators were attached to the uppermost cylinder: the three representing Vastus Lateralis (VL), Vastus Medialis (VM), Vastus Intermedius (VI), attached to the femur, and the one representing the Rectus Femoris (RF) attached on the pelvis. The hamstrings were composed of Semitendinosus (ST), Semimembranosus (SM), Biceps Femoris long head (BF-lh) and Biceps Femoris short head (BF-sh). These two last muscles converged in a single insertion point on the head of the fibula. Their origin was respectively on the ischial tuberosity (BF-lh) and the lateral side of the femoral shaft (BF-sh). The SM and ST originated in a narrow area of the ischial tuberosity and were connected to the medial

side of the proximal tibia. A via point was defined for the SM to reproduce the deviation of the action line of this muscle occurring when the knee joint approaches the maximum extension. The Triceps surae muscle was considered composed of Gastrocnemius Lateralis (GaL), Gastrocnemius Medialis (GaM) and Soleus. Only the two gastrocnemii affect the knee joint, acting as knee flexors. They were represented by having an insertion point on the rear side of the lateral and medial condyles respectively and converging to a common attachment point located on the rear of the foot (the Achille's tendon, on the heel).

#### *2.2. The Walking Model*

To assess the performance of the knee model during walking, a driving model was implemented. It was composed of trunk, pelvis, thighs, shanks and feet (Figure 3) and was animated by kinematic data reproducing a natural walking. The femur of the knee model was attached to the thigh segment of the walking model, so that it underwent the same accelerations/decelerations. A geometric solid representing the foot was attached to the distal tibia.

**Figure 3.** The walking model. It was constituted by the trunk (not represented), the pelvis, the thighs, the shanks and the feet. The left limb is represented with bones and muscles. Sample frames from the stance phase are reported above, sample frames from the swing phase in the second row. HS: heel strike; Cnt-TO: contralateral toe off; Cnt-HS: contralateral heel strike; TO: toe off. All joints were controlled by rotation motors: 3 degrees of freedom (d.o.f) between trunk and pelvis, 3 d.o.f. between pelvis and thigh, 2 d.o.f. between thigh and shank, 2 d.o.f. between shank and foot.

All segments of the walking model were considered as rigid bodies and were linked together by rotational joints.

The trunk was taken as the reference segment. The pelvis was attached to the trunk by a series of rotation actuator (motors) that controlled the three degrees of freedom (d.o.f.) of the pelvis: anteversion/retroversion, frontal plane tilt and horizontal rotation. Similar constraints controlled the 3 d.o.f. of the thigh about pelvis: hip flexion/extension; adduction/abduction; internal/external rotation; the 2 d.o.f. between shank and thigh, knee flexion/extension and internal/external rotation; and the 2 d.o.f. between foot and shank, ankle plantar/dorsiflexion and pronation/supination.

The walking model was animated by data from our repository (Movement Biomechanics and Motor Control Lab, Politecnico di Milano, the SAFLo (Servizio di Analisi della Funzionalità Locomotoria) protocol [37–39]). These data referred to the average of 14 strides from 5 normal subjects, males, age between 24 and 36 years, walking barefoot at their natural velocity. The adopted acquisition protocol [38,39] provided the kinematic

data corresponding to the above-defined degrees of freedom. In addition, the average ground reaction force and its point of application were available. All variables were normalized in time and referred to an ideal cycle duration of 1 s. Masses of body segments, corresponding to an ideal subject 1.72 m tall and 70 kg of body mass, were obtained from an anthropometric table [40]. The inverse dynamics problem solution provided us with the knee joint moments corresponding to this ideal subject. When the knee model was attached to the walking model, care was taken to make the femoral head to coincide with the center of the acetabulum in the pelvis. During the animation of the walking model, the flexion/extension was imposed on the knee model while the remaining five degrees of freedom were let unconstrained. This was performed by implementing a special device that we called "the Grood&Suntay mechanism". " "

#### *2.3. The Grood&Suntay Mechanism*

The principle was based on the convention adopted by these authors [41] to uniquely define the functional axes and the six degrees of freedom of the knee. The mechanism is described in Figure 4. The uppermost shaft was connected to the femur. A revolute joint, positioned in the center of the knee, allows for a flexion/extension (F/E) movement of a plate (red in the figure). A second revolute joint, with an axis perpendicular to the previous one, allows for the adduction/abduction (Ad/Abd) movement of a second plate (green) with respect to the red one.

**Figure 4.** Left: the so-called Grood&Suntay mechanism (G&S) used to measure the six d.o.f. of the knee and to impose the flexion/extension movement (F/E). A generic position is represented, with the tibia translated in relation to the center of rotation. Right panel, up: the forces applied to the upper tibia deriving from the ground reactions; middle: the corresponding transfer moments (values are positive for flexion, external rotation, adduction moments); bottom: the flexion/extension moment applied by the G&S mechanism in the absence of muscle forces (negative means extensor).

A third revolute joint has a rotational axis perpendicular to the previous one and connects the green plate to the lowermost shaft, which will be aligned with the tibia longitudinal axis. This revolute joint will permit the internal/external rotation (I/E rot) of the shank. The three rotational axes converge in a single point that can be considered the center of the knee joint, and thus, the mechanism corresponds to a spherical joint,

of which the three rotational components are defined. However, the tibia should not be constrained by a spherical joint (we want to analyze also the translational movements), and so the connection between the lowermost shaft and the tibia was designed so that the shaft can translate in the three directions (X, Y, Z in Figure 4) by keeping the same orientation with respect to the tibial axis. A rigid plate (grey in the figure) was rigidly attached to the proximal tibia, as to provide a reference for the translations, and the "parallelism" constraint was applied between the lowermost shaft and this plate. In this way, the three components of the relative translation between tibia and femur can be measured. They are the anterior/posterior (Z), medial/lateral (X), and distal/proximal (Y) components. Besides providing the six free components of the tibia-femur motion, the mechanism can also be used to drive the knee joint kinematics. In our case, we inputted the time course of the flexion/extension angle to the corresponding revolute joint, and the flexion/extension movement was obtained without any constraint for the remaining five degrees of freedom.

#### *2.4. Loading of the Knee Model*

The external loads applied to the knee model were the ground reaction forces measured during walking and the muscle forces. The inertia components (moments and forces) associated to the movement of shank and foot were computed by applying an inverse dynamics analysis. Concerning the ground reaction forces (anterior/posterior, medial/lateral, vertical components), we deemed necessary to transfer them to the center of the knee joint because when they were applied to the foot, which was attached to the knee model, the moment they produced at the knee was affected by the position of the foot, which in turn was dependent on the knee moment, and thus, the system became unstable. When we moved the ground reaction force to the upper extremity of the tibia (the weights of foot and shank were subtracted to the vertical component), the corresponding transfer moment was applied to the tibia. As reported in Figure 4, the transfer moment is null during the swing phase because it is the moment that originally the ground reaction forces produced at the knee joint when they were applied in foot/ground contact area. Instead, the moment required to produce the knee flexion/extension includes the inertia components as well. In the absence of muscle forces, this moment was produced by the G&S mechanism. Its time course is reported in the right lowermost panel of Figure 4. As to the muscle forces, an optimization procedure was implemented to obtain an estimate of the muscle forces required to sustain the external knee joint moment. The criterion was the minimization of the maximum force (Min/Max criterion as suggested by [42]) in relation to the physiological cross-sectional area of the muscle (PCSA). For this computation, the muscles considered were not only those acting at the knee but also Iliacus, Psoas, Glutei (maximum, minimum and medium), Adductors (magnus, longus and brevis), Tensor fasciae latae, Soleus, Tibialis anterior and posterior and Peroneus, for a total of 23 muscles. The PCSA values for all the muscles were obtained from tables reported in [43]. The optimization was implemented with the following constraints:

$$F\_l \ge 0 \; ; \quad \sum\_i F\_l b\_{i\bar{j}} = M\_{\bar{j}} \tag{3}$$

where *F<sup>i</sup>* is the force of the muscle *i*; *Mij* is the moment produced by muscle *i* around the functional axis *j*; *bij* is the lever arm of muscle *i* in relation to axis *j*.

The lever arms of the muscles could not be obtained by simple geometric considerations, because the lines of action of the force actuators representing the muscles changed along the movement, the tendons of some muscles wrapped over the bone surfaces, and these surfaces were not always aligned with the rotation axes. The lever arms were thus obtained by analyzing the moment produced by each muscle when it was activated with a predefined muscle force. For each muscle, a simulation with no gravity, no ground reaction forces and no mass of the segments was run. All muscle forces were set to zero except for the muscle considered. A constant force of 100 N was imposed on this muscle, and the resulting moments (flexion/extension, adduction/abduction and internal/external

rotation) were analyzed. At each joint (hip, knee and ankle), the lever arm corresponding to each rotational degree of freedom was computed as *bij* = *Mij* . = 

*Fi* Figure 5 reports the estimated muscle forces referring to the muscles of interest for the knee joint. 

 ≥ 0 ; ∑ = 

**Figure 5.** Muscle forces applied to the model by the muscle actuators. Dashed lines represent the sum of the muscle forces belonging to the respective muscle group: Quadriceps (Quadr), Hamstrings (Hamst), Gastrocnemius (Gastr).

#### *2.5. Simulation Conditions*

The walking model received the kinematic variables as input. The knee model was attached to the walking model through a rigid connection between femur and thigh. The only flexion/extension was imposed on the knee, while the remaining degrees of freedom were determined by forward dynamics computation. The Kutta–Merson integration algorithm was adopted with rendering frequency of 100 frames/s and integration step ranging from 0.010 to 0.001 s; configuration tolerance was about 0.001 mm for displacements and 0.05 ◦ for rotations. To understand the role of the knee morphology and kinematics over the tension of the ligaments, a reference simulation was run without external loads and muscle forces. Then the ground reaction forces, the moments and the muscle forces were applied to the knee model, and the changes in ligament tensions were analyzed. In a second series of simulations, the forces of one of the three muscle groups were sequentially set at zero: the quadriceps muscles (NoQuad), the hamstrings muscles (NoHam) and the gastrocnemius muscles (NoGas). This allowed us to infer about the effect of each muscle group on the tension of the ligaments. The changes produced by the elimination of each muscle group were quantified by measuring the average difference of the ligaments' tensions between the "All muscles" condition and the condition in which a single muscle group was removed. The analysis was carried out in the phases in which each muscle group was dominating: the first contraction of quadriceps occurred from 0% to 30% of the stride cycle, and this phase was named the Loading Phase; from 31% to 80% the second contraction of quadriceps occurred and was named the Knee Flexion Phase; from 81% to 100% the dominating muscles were the hamstrings, and this phase was named Knee Breaking Phase; from 20% to 60% of the stride cycle the gastrocnemius muscles contraction occurred and the phase was named the Support Phase.

#### *2.6. Model Validation*

A qualitative assessment of the performance of the model along the whole stride cycle was accomplished by looking at three major determinants of the knee kinematics: (a) the forward displacement of the tibia during knee flexion, (b) the presence of the screw home mechanism (external rotation of the tibia when approaching the full extension and internal rotation associated to the initial knee flexion), (c) the time course of the tibio-femoral contact force that must exhibit two peaks during the stance phase (the second higher than the first, [15,44]) and never become zero during the swing phase, in order to guarantee the joint stability. The quantitative assessment was made by comparing our data with those provided as a reference for the Grand Challenge competition [17], which are available on line at the link: https://simtk.org/projects/kneeloads (accessed date 15 January 2021). The data about normal walking from the 2012 competition were downloaded and analyzed.

#### **3. Results**

The comparison between our data and the reference data showed that our tibiofemoral contact force in the "All muscles" condition was larger than the reference one (see Figure 6, up). The main features of the two curves, however, were similar: two peaks in the stance phase, the second higher than the first one, and a third peak at the end of the swing phase. The root mean square (RMS) of the difference was 429 N. Downscaling the contact force by a factor 0.7 reduced the RMS to a minimum of 318 N. The lower panels of Figure 6 report the knee joint kinematics as measured by our G&S mechanism in the two loading conditions: "All muscles", in which all external loads and muscle forces were applied to the model, and "No muscles", in which the knee joint was moved without muscle forces (ground reactions and moments were also removed). In this condition, the rotations and the displacements were determined by the only contact surfaces and by the ligaments' tensions. Their amplitudes were in the range described in the literature, both for the loaded and the unloaded conditions [45,46]. The amplitude of the internal/external rotations, anterior/posterior displacements and medial/lateral displacements were slightly increased when loads, and muscle forces were applied, although the general time course was not very different. No appreciable difference appeared concerning the adduction/abduction and the proximal/distal movements. Some oscillations were observed all along the stride cycle, more pronounced in the "No muscles" condition. They could have been triggered by the residual roughness of the contact surfaces. Nevertheless, the physiological screw-home mechanism (knee external rotation associated with knee extension and vice versa) was evident either in the unloaded and in the loaded condition. As to the displacements, a forward displacement of the tibia was observed during both knee flexion phases (the first during stance and the second associated with the swing phase). The medial/lateral displacement was also in relation to the two flexion phases of the knee. The forces supported by the cruciate ligaments and the collateral ligaments are reported in Figure 7. As to the cruciate ligaments, it appears that the muscle activity produces a much higher tension (continuous lines) than in the "No muscles" condition (dashed lines). This does not happen for the collateral ligaments, where the tensions in the two conditions are not very different. The ACL is loaded in two phases along the stride cycle, the first in correspondence with the load acceptance phase and the second in correspondence with the initiation of the unloading phase, when the knee starts flexing. The PCL presents a mild tension in the stance phase, in correspondence with the knee extension movement, and considerably high tension in the

second half of the swing phase, in correspondence with the large knee extension movement. For both cruciate ligaments, there seems to be no specific engagement angle but rather a strong dependence on loading conditions.

displacements of the tibia in relation to the femur as measured by our "Grood&Suntay" reported: "All muscles" refers to the application of muscular forces simulating a normal loading condition, and "No muscles" refers to the **Figure 6.** Up: comparison of tibio-femoral contact forces between our result (**blue**) and reference data (**black**), taken from the Grand Challenge competition (see text). The orange curve represents our contact force after downscaling by a factor of 0.7. Mid and lower panels: rotations and displacements of the tibia in relation to the femur as measured by our "Grood&Suntay" mechanism (G&S). Flexion/extension was imposed, all other movements resulted from the dynamic equilibrium of forces and by the geometry of the contact surfaces. Two conditions are reported: "All muscles" refers to the application of muscular forces simulating a normal loading condition, and "No muscles\*" refers to the absence of muscle forces and ground reactions.

" " " **–** " " " **Figure 7.** Tension of the cruciate ligaments (**up**) and collateral ligaments (**down**) along the stride cycle. Continuous lines refer to the "All muscles" condition, dashed lines refer to "No ground reaction forces–no muscles forces" condition. For this condition, the muscle labels are evidenced by an "\*".

**–** The two collateral ligaments (Figure 7, lower panel) present a baseline tension around 20–25 N. The MCL is recruited in three phases: early stance, late stance, early swing. In the second half of the swing, its tension drops below the baseline. The LCL instead is loaded during the whole swing phase. All components of MCL and LCL (MCL-a, MCL-i, MCL-p, LCL-1, LCL-2, LCL-3, not represented in the figure) had similar behavior.

The tension of the deep components of MCL and the capsular ligaments are reported in Figure 8. The MCL-dp was loaded during the stance phase, with two peaks in correspondence with maximum knee extension, while the anterior component (MCL-da) was loaded in the swing phase.

" "

" " " **–** " " " **Figure 8.** Tension of the deep medial collateral ligament (**a**) and capsule ligaments along the stride cycle (**b**). Continuous lines refer to the "All muscles" condition; dashed lines refer to "No ground reaction forces–no muscles forces" condition. For this condition, the muscle labels are evidenced by an "\*".

" "

" "

The posterior components of the capsule (Cap-Post-L and Caps-Post-M) were also recruited during the stance phase: the Cap-Post-L more continuously and the Cap-Post-M with two peaks at the early stance and late stance, respectively, and a zero tension in the middle. The tension of the anterolateral capsule (Cap-Ant-L) occupied the whole stance phase and half of the swing phase. All components of the capsule exhibited considerably reduced tension in the "No muscles" condition. " " " **–** " " "

When the simulations were carried out by sequentially removing the quadriceps or the hamstrings, the obtained results were the ones reported in Figure 9. Compared to the "All muscles" condition, the ACL exhibited a drastic reduction of tension when the quadriceps was removed. Only a short phase of loading remained at the beginning of the stride cycle. The removal of hamstrings instead did not affect the ACL tension, except at the end of the swing phase, where the tension appeared enhanced. Removing the quadriceps had an opposite effect on the PCL tension: The phase of recruitment during the stance phase was anticipated and lasted for the whole stance phase, with the appearance of two peaks considerably higher than in "All muscles" condition. The PCL recruitment in the second half of the swing phase was almost unaffected by the quadriceps removal. At contrary, the "NoHam" condition did not affect the PCL tension in the stance phase but drastically reduced the PCL tension in the last swing phase. " " " " " "

**Figure 9.** Changes of ligaments tensions obtained by removing the quadriceps group (NoQuad) or the hamstrings group (NoHam). The black curve represents the reference condition in which all muscle forces were applied in the model.

Concerning the collateral ligaments (Figure 9, lower panel), the removal of quadriceps muscle produced a generalized reduction of tension in both LCL and MCL. The appearance of a short-lasting recruitment of LCL ligament at around 0.2 s is probably related to an abnormal position of the two contact surfaces that was caused by the unbalanced muscular condition. The removal of the hamstrings produced no relevant changes except at late swing phase, where the LCL tension appeared reduced and the MCL tension appeared increased.

The elimination of quadriceps had a mild effect on the deep components of the MCL (see Figure 10): only in a small portion of the recruitment phase of MCL-dp the tension was appreciably reduced. The quadriceps removal had instead a considerable effect in reducing

the force in all the components of the capsule: Cap-Ant-L, Cap-Post-L, and Cap-Post-M. No effect was produced by the elimination of the hamstrings, except for a short-lasting increase of the tension in the Cap-Ant-L at late swing phase.

**Figure 10.** Changes of ligaments tensions obtained by removing separately the quadriceps group (NoQuad), the hamstrings group (NoHam) or the gastrocnemius group (NoGas). The black curve represents the reference condition in which all muscle forces were applied in the model.

" " The effects of the elimination of the gastrocnemius force were not included in the previous figures for the sake of clearness. Just as an example they are reported in Figure 10 with reference to the three capsular components. However, the tension of the ligaments changed also in the PCL and LCL components as it will be shown later.

The quantification of the effects of muscles removal is reported in Table 2. The difference of ligaments tension between "All muscles" condition and the absence of one muscle group was averaged by keeping separately the phases in which the difference was negative (tension reduction) from the ones in which it was positive (tension increased). The opposite behavior of ACL and PCL is evident: the elimination of quadriceps produced a reduction of tension in the ACL and an increase in the PCL, while the elimination of hamstrings produced the increase of ACL tension and a decrease in the PCL. The NoGas condition produced an unclear situation in the ACL: positive and negative variations were almost equivalent, but a net prevalence of tension increase in the PCL was observed. Concerning the other ligaments, the NoQuad condition produced in both load acceptance and knee flexion phases a prevalence of tension reduction in all MCL and capsular ligaments. Interestingly, in the same NoQuad condition, the LCL tension increased in the load acceptance phase and decreased in the knee flexion phase. LCL tension decreased also in the NoHam condition and increased when the gastrocnemius muscles were removed. Other remarkable variations were the increase of MCL and anterior capsule tension when the hamstrings were removed, and the decrease of capsular ligaments tension with gastrocnemius muscles removal. No clear effect could be observed on the deep components of MCL, in which the negative and positive changes were almost equivalent in all muscle removal conditions.


**Table 2.** The average difference of ligaments tensions (N) between "All muscles" condition and removal of one muscle group. The averages were computed separately for phases in which the difference was negative and phases in which the difference was positive. The light blue and the pink evidenced cells represent the predominant effect (negative or positive variations, respectively) resulting from the removal of the corresponding muscle group in the respective phase.

#### **4. Discussion**

The objective of our work was to investigate how the muscle contractions occurring during walking affect the knee ligaments tension. The problem is of interest because of its implications with surgery planning (ligaments reconstruction and knee joint arthroplasty) and rehabilitation (recovery of muscle force, compensation strategies and orthoses). Although the role of these ligaments is well established in general terms [34], only a few publications deal with knee ligaments behavior in living subjects and dynamic conditions, [22,23,35,47–49] and only one, to our knowledge, presents an estimation of the tension of the ligaments along a stride cycle, based on a musculoskeletal model [27]. The reason could be related to the complex structure of the knee joint and to the difficulty to predict the real loading conditions. Our musculoskeletal model represents a compromise between the need to represent that structure in deep detail and the need to realize a flexible tool that can simulate different functioning conditions. The performance of the model was relatively good. Even if the tibio-femoral contact force predicted by our model was larger than the force measured by an instrumented prosthesis [15,16,50], the main features were well reproduced. About the amplitude, it must be considered that our input data are from a different source, and probably, the anthropometric parameters and the walking velocity were different. In fact, after downscaling the absolute values, the RMS of the difference with respect to reference data was 318 N. Let us consider that the winners of the Grand Challenge competition [17] estimated the tibio-femoral contact forces with a RMS in the range from 229 N to 312 N for medial contact force and from 238 N to 326 N for lateral contact force. The main limits of our model are the lack of an interface between bone surfaces representing the cartilage and the menisci, the simple representation of the patellofemoral joint and the difficulty to adapt the model to different sizes, mainly because the repositioning of the ligaments is a very delicate procedure. Therefore, at present, it can be considered a generic model that hardly could become subject specific as it would be desirable. As to the muscle forces that are inputted to the model, different optimization approaches could have been used. Most of them are based on hypotheses about the motor control

system [2,51–53], and their results are clearly dependent on the cost function adopted and on the many parameters that have to be assumed (see [54] for an overview of the methods for muscle force prediction). For this reason, there are inconsistencies among the different authors [31,55–57]. In our case, we were not interested in the investigation of motor control strategy, so our approach was a simple one, though relatively efficient [42], that just looks for the minimization of the maximum muscle force referred to the physiological crosssectional area of the muscle. The results are consistent with the known electromyographic patterns [58], as for most of the published works [31,51,52,56,57]. We suppose that even if the estimation of the muscle forces had been slightly different, the effect of these forces on the tension of the ligaments would not have changed substantially. Concerning the biomechanical structure, while designing our model, we realized how difficult was to arrange the ligaments in the model to avoid the excessive resistance to movement and, at the opposite, to prevent joint dislocation. Beyond the generic anatomical description of the attachment areas [34], the specific data reported in the literature [6,22,27,59] did not help very much: they showed that the attachment points can differ among subjects, and a great variability exists about the mechanical properties. The matching between attachment points, ligaments properties and bone surface geometry is very critical. For example, by using the force/deformation characteristics and stiffness values provided in [26], the anterior bundle of the anterior cruciate ligament (ACL-a) reaches a force of 45 N for a deformation ε = 0.06, which correspond to a lengthening of just 2.1 mm (rest length L<sup>0</sup> = 35.2 mm). A stiffer ligament, for example the deep posterior bundle of the medial collateral ligament (MCL-dp), goes from 0 N to 135 N for a length increase of just 2.4 mm and achieves 217 N for elongation of 3 mm. Therefore, the acceptable range of length changes is very small. Furthermore, the force distribution between all ligaments affects the knee joint kinematics. In our model, the only flexion/extension movement was imposed, and the remaining 5 d.o.f. were completely free. Obtaining the results reported in Figure 6 was an important achievement in our point of view. In fact, the rotations and the displacements are within a reasonable range and the screw home mechanism is reproduced both in the unloaded condition and also when the external forces and muscle forces were applied. About the muscle forces, our estimation takes into account the muscles arrangements, their lever arms in relation to the different functional axes and the physiological cross-sectional areas. A total of 23 muscles were considered because some of the knee muscles are double-joint muscles, acting also at the hip (hamstrings) and the ankle (gastrocnemius muscles), and so the optimization procedure has to take into account the whole limb. A comparison of our result to a reference pattern of muscle force generation is difficult, because, as mentioned above, there are differences between the data provided by the different authors. However, the phases of activity and the relative amplitudes obtained for our muscles (see Figure 5) are consistent with their role, and can easily be understood: the first phase of the activity of the vasti group (from 0 to 30% of the stride cycle) corresponds to the load acceptance phase, in which the quadriceps is required to control the knee flexion after the impact with the ground; the second phase of activity, starting in the second half of the stance phase, is required to sustain the knee flexion while the ground reaction is still present. The activity of the vasti is extended to the first part of the swing phase, since an extensor moment is required to prevent excessive knee flexion while the limb is thrown forward. Interestingly the rectus femoris (RF) becomes active in this phase, since its action is required at the hip joint to produce the forward acceleration of the thigh, and thus, it contributes to the extensor moment together with the vasti muscles.

At the end of the swing phase, the hamstrings are the muscles required to reduce the forward velocity of the limb, and this is achieved by producing an extensor moment at the hip and a simultaneous knee flexion moment. The gastrocnemius muscles force dominates most of the stance phase since the major role of this muscle is to produce the plantarflexion moment at the ankle, synergistically with the soleus.

When these muscle forces and the ground reaction forces were applied to the knee, the kinematics changed with respect to the unloaded condition, but still preserved its basic

features, as represented in Figure 6. A notable difference was observed in the backward displacement of the femur, which was interrupted at about 75% of the stride cycle in the loaded condition. This is the effect of the activity of the hamstrings that brake the forward movement of the shank and so produce a rear displacement of the proximal tibia in relation to the femur. The hamstrings activity is also responsible for the increase of tibio-femoral contact force that is observed at the end of the swing phase, both in our data and in the reference data (Figure 6). The adduction/abduction rotation and the distal/proximal displacement were relatively unaffected by the muscle contractions. The changes in the knee kinematics had as a counterpart a different loading of the ligaments. Considerable was the effect on the cruciate ligaments. The total force of the ACL increased, in the first phase, corresponding to early stance phase, from 20 N to about 200 N, as well as in the second phase, corresponding to late stance-early swing the force changed from 0 N to about 100 N. The maximum PCL force was about 150 N in the unloaded condition and occurred at around the 80% of the stride cycle. In the loaded condition, the PCL tension increased up to 670 N with a peak at approximately 90% of the stride cycle. The relation between the increase of tension in this phase and the presence of hamstrings force, which produce a rear displacement of the tibia is evident. It is interesting to note that in the unloaded condition the ACL was recruited at around knee maximum extension (about 15◦ of knee engagement angle) and the PCL had a maximum tension in correspondence with maximum knee flexion, while in the loaded condition the engagement of the cruciate ligaments was not closely related to the knee angle, but it was strictly dependent on the forces applied. The different components of the collateral ligaments were less affected by the application of the muscular forces, and this is consistent with their longitudinal arrangement in relation to the long axis of the limb. Actually, the distal/proximal displacement and the adduction/abduction movements could only have been changed if the femoral condyles and tibia surfaces detached from each other, in one or both knee joint compartments, which did not occur. Therefore, the small increase of the strain (and consequently the change in ligament tension) was in this case mainly due to the increase of internal/external rotation. Both anterolateral (CAP-Ant-L) and posterior capsule (Cap-Post-L and Cap-Post-M) were affected by muscle contractions and almost doubled their tension. This is evidently associated to the increased internal/external rotation and anterior slide of tibia produced by the quadriceps contraction in the same phase of the walking cycle.

The effect of the muscle contractions on the tension of the ligaments became more evident when the single muscle groups were removed in our simulations. The most affected ligaments were again the cruciate ligaments. ACL tension dropped dramatically when quadriceps was eliminated, thus revealing that quadriceps was responsible for ACL loading. PCL force instead increased with quadriceps force removal, which indicates that the remaining forces applied to the knee produced a tendency for the proximal tibia to slide backwards. Removal of the hamstrings produced an increase of the ACL tension at the end of the swing phase and a reduction of the tension in the PCL. Removal of the gastrocnemius muscles had a not well-defined effect on the ACL (the negative and positive differences were almost equivalent) and a positive effect on the PCL and LCL. All the effects of muscle removals in the other ligaments are well depicted in Table 1. These results could only be compared to similar results referring to natural walking [27,60]. The data reported in the first of these papers seem difficult to interpret though, in that the ACL appears constantly loaded for the whole stance phase, part of the swing phase and again at late swing; the contribution of LCL and MCL is mild, and the PCL appears unloaded in the stance phase and only lightly loaded at the end of the swing phase. Unfortunately, the second paper does not present the time course of the ligaments tension nor the muscle forces inputted to the model. Our results thus can only be considered in relation to the general clinical knowledge eventually supported by a biomechanical analysis [8,22,49,61,62]. All these authors argue that ACL tension strongly depends on quadriceps contraction and that hamstrings activity can help to reduce the ACL load in subjects affected by ACL deficiency. The fact that our results are perfectly in line with those observations gives us some confidence about the reliability of our model and its ability to effectively represent the knee joint biomechanics.

#### **5. Conclusions**

The analysis of the tension of the ligaments in the knee joint is a complex problem that has been faced through different approaches, ranging from studies on cadaver specimens, to direct measurement by implanted sensors to biomechanical modelling. Few studies deal with natural walking, and the results are still questionable. Our musculoskeletal model study provides a contribution to understanding the role of the different muscles in determining the ligaments loads. Although the method would require further validation, the main consolidated phenomena about the knee joint mechanism, like the screw home mechanism, the posterior displacement of the femur during flexion, the loading/unloading effects on the cruciate ligaments produced by the contraction of quadriceps and hamstrings, were confirmed in our simulations. Hence, the model appears to be an effective tool for further investigating the biomechanics of knee joint under dynamic conditions.

**Author Contributions:** Conceptualization and methodology, C.A.F.; software implementation and simulations C.A.F. and L.D.; data curation C.A.F. and L.D.; original draft preparation C.A.F.; review and editing C.A.F. and L.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** This study did not require ethical approval, as no experimental data was acquired specifically for this study.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The contribution of Maddalena Grossi to the estimation of muscle forces through the optimization approach is warmly acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Biomechanical Consequences of Tibial Insert Thickness after Total Knee Arthroplasty: A Musculoskeletal Simulation Study**

**Periklis Tzanetis 1,\* , Marco A. Marra <sup>1</sup> , René Fluit <sup>1</sup> , Bart Koopman <sup>1</sup> and Nico Verdonschot 1,2**


**Featured Application: The proposed methodology provides orthopedic surgeons with quantitative data on the sensitivity of muscle, ligament, and joint compressive forces to tibial polyethylene insert thickness variations and can be further developed to optimize navigated or roboticassisted total knee arthroplasty.**

**Abstract:** The thickness of the tibial polyethylene (PE) insert is a critical parameter to ensure optimal soft-tissue balancing in the intraoperative decision-making procedure of total knee arthroplasty (TKA). However, there is a paucity of information about the kinetic response to PE insert thickness variations in the tibiofemoral (TF) joint, and subsequently, the secondary effects on the patellofemoral (PF) biomechanics. Therefore, the purpose of this study was to investigate the influence of varying PE insert thickness on the ligament and TF compressive forces, as well as on the PF forces and kinematics, after a cruciate-retaining TKA. A previous patient-specific musculoskeletal model of TKA was adapted to simulate a chair-rising motion in which PE insert thickness was varied with 2 mm increments or decrements compared to the reference case (9 mm), from 5 mm up to 13 mm. Greater PE insert thickness resulted in higher ligament forces and concurrently increased the TF compressive force by 21% (13 mm), but slightly unloaded the PF joint with 7% (13 mm) while shifting the patella distally in the trochlear groove, compared to the reference case. Thinner PE inserts showed an opposite trend. Our findings suggest that the optimal PE insert thickness selection is a trade-off between the kinetic outcomes of the TF and PF joints.

**Keywords:** musculoskeletal model; polyethylene thickness; total knee arthroplasty; ligament forces; compressive forces; patellar kinematics

#### **1. Introduction**

Total knee arthroplasty (TKA) is an effective surgical intervention for end-stage knee osteoarthritis in which the diseased articulating surfaces of the knee joint are replaced with artificial components to relieve pain and restore the normal knee function [1]. A successful surgical outcome depends, amongst others, on the design [2], size [3], and alignment of the prosthetic components [4]. A common mode of early surgical failure includes mechanical instability that may result from aseptic loosening or polyethylene (PE) wear [5]. These failure mechanisms have been associated with the inadequate thickness of the tibial PE insert [6]. The selection of an optimum PE insert thickness in the intraoperative procedure is important to achieve soft-tissue balance and restore the physiological function in the tibiofemoral (TF) joint, thus lessening the risk of implant-related complications [7]. Previous studies suggest that a minimum thickness of 6 to 8 mm is required to minimize the contact stresses within the tibial insert surface [8,9], applying the conventional 2 mm increments; these increments allow surgeons to identify the effect of PE insert thickness

**Citation:** Tzanetis, P.; Marra, M.A.; Fluit, R.; Koopman, B.; Verdonschot, N. Biomechanical Consequences of Tibial Insert Thickness after Total Knee Arthroplasty: A Musculoskeletal Simulation Study. *Appl. Sci.* **2021**, *11*, 2423. https:// doi.org/10.3390/app11052423

Academic Editor: Carlo Albino Frigo

Received: 4 February 2021 Accepted: 4 March 2021 Published: 9 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

variations on the soft-tissue tension and intraoperative kinematics based on knee laxity trials [10,11]. Following this recommendation, some authors argue that the insertion of thicker PE components is likely to prevent accelerated wear and early surgical failure [12]. Nevertheless, the debate continues about the target PE insert thickness, with contradictory findings from earlier studies on the topic. Greco et al. [13] demonstrated that implanted knees with a PE insert thickness of more than 14 mm performed similarly in the postoperative clinical evaluation, yielding no surgical failure due to aseptic loosening or tibial component instability as compared to those with a thickness under 14 mm. Contrarily, Berend et al. [14] labeled tibial PE inserts thicker than 14 mm as potential contributors to posterolateral rotatory instability, leading to higher failure rates at mid- to long-term follow-up. One of the main limitations of these studies, however, was the failure to control for the variability in the depth of tibial resection, which might provide an inaccurate estimation of the actual effect of PE insert thickness variation on the surgical outcomes.

From a biomechanical perspective, using thicker tibial PE inserts than the minimum recommended threshold (6–8 mm), while keeping the resection depth unchanged, can increase tension in the soft-tissue envelope; this can subsequently result in higher contact forces in the distal femur and tibial insert interface. The increased contact forces are also transferred to the fixation site, and this could, therefore, promote aseptic loosening of the prosthesis. A comprehensive evaluation of the kinetics after TKA is crucial to prevent higher than physiological mechanical stress and strain on the ligamentous structures as well as high load levels at the interface between bones and prosthetic components. This knowledge could lead to better surgical choices and increased implant longevity. Since the thickness of the tibial PE insert can be modified intraoperatively, the surgeon can customize the surgical plan by choosing a PE insert thickness that best replicates the natural knee kinematics and kinetics. Computational models of the musculoskeletal system represent a valuable tool in this intraoperative decision-making procedure. Musculoskeletal models can be personalized to patient-specific anatomical features and predict variations in the knee joint loads and kinematic patterns due to changes in implant design properties or surgical choices, such as the thickness of the tibial PE insert. These predictive capabilities of musculoskeletal models provide quantitative patient-specific information that can aid surgeons to formulate an optimal treatment plan tailored to individual patients [15].

Earlier studies typically focus on the primary effects of PE insert thickness on the TF joint biomechanics. For instance, researchers have recently examined the effect of varying PE insert thickness on the TF kinematics after TKA, using a dynamic knee simulator of laxity tests [16], or by performing an intraoperative assessment using a computer-assisted surgery navigation system [7]. However, to the best of our knowledge, current scientific literature lacks a proper investigation of the relationship between PE insert thickness and TF kinetics, as well as of the effects on the patellofemoral (PF) joint where the secondary effects may influence the biomechanical parameters, including forces and kinematics, in a clinically relevant way. Therefore, the aim of this study was to first investigate the influence of PE insert thickness, in isolation, on the TF joint in terms of ligament and compressive forces, and secondly, to assess to what extent PE insert thickness variation affects the biomechanical parameters on the PF joint. For this, we used a patient-specific musculoskeletal model with a cruciate-retaining prosthesis. We hypothesized that thicker tibial PE inserts would result in higher ligament forces and, thus, increased TF compressive force, but the effect would remain neutral with regard to the PF articulation.

#### **2. Materials and Methods**

A previously validated musculoskeletal model of a TKA patient was used in this study [17]. In brief, the model was built upon the Twente Lower Extremity Model 2.0 (TLEM 2.0) template [18] and integrated patient-specific data available from the fifth "Grand Challenge Competition to Predict In Vivo Knee Loads" dataset [19]. This lowerextremity model of the implanted side comprises 11 degrees of freedom in the TF and PF joints. Flexion and extension angle in the TF joint is motion capture driven, whereas

the kinematics of the remaining degrees of freedom are solved quasi-statically, using the force-dependent kinematics (FDK) approach, as proposed by Andersen et al. [20]. This methodology enables the concurrent estimation of muscle, ligament, and joint contact forces.

For this study, we introduced some changes with respect to the existing patientspecific model, using the AnyBody Modeling System v. 7.3.1 (AnyBody Technology A/S, Aalborg, Denmark); software for modeling and simulation analysis of the musculoskeletal apparatus [21]. The patellar tendon was modeled as three non-linear spring elements instead of a single rigid element connecting the patella and tibial tuberosity. Stiffness of the patellar tendon was determined by multiplying the median cross-sectional area from the patient-specific computed tomography (CT) with a modulus of elasticity reported in the literature [22]; likewise, stiffness and reference strain of the TF and PF ligament bundles were adopted from the existing musculoskeletal model configuration. In addition, we configured the model to incorporate a size 6 Triathlon (Stryker, Kalamazoo, MI, USA) cruciate-retaining implant (Figure 1), ensuring that the femoral anteroposterior anatomical dimension prevents anterior notching of the femoral cortex. This implant was chosen as it has shown good functional performance and ligamentous stability [23]. The femoral and tibial component geometry was positioned in such a way as to match tangentially the patient-specific distal femoral and proximal tibial bone cuts, applying the well-established mechanical alignment technique [24]. The patella was not resurfaced. Since the cartilage underside of the patella could not be determined from the patient-specific dataset, we generated a 3 mm offset on the dorsal facet of the patellar bone according to the mean patellar cartilage thickness estimation of Cohen et al. [25]. In this model, PE insert thickness was varied with 2 mm increments or decrements (−4 mm, −2 mm, +2 mm, +4 mm) compared to the reference case (9 mm), testing five thickness cases in total: 5 mm, 7 mm, 9 mm, 11 mm, and 13 mm (Figure 2). The effect of PE insert thickness was evaluated by simulating a chair-rising activity based on motion capture data available as part of the open access grand challenge competition dataset.

**Figure 1.** Illustration of the patient-specific musculoskeletal force-dependent kinematics (FDK) knee model implanted with a Triathlon cruciate-retaining prosthesis. This full-body model includes the head, trunk, pelvis, and the lower extremity of the implanted side (thigh, shank, patella, talus, foot) (left). A close-up anterior view of the knee joint with the medial collateral ligament (MCL), lateral collateral ligament (LCL), posterior cruciate ligament (PCL), medial patellofemoral ligament (MPFL), lateral patellofemoral ligament (LPFL), and patellar tendon (PT) (right).

**Figure 2.** Illustration of three custom postoperative cases simulated in this study. The tibial polyethylene (PE) insert thickness cases of 9 mm (reference), 11 mm (+2 mm), and 13 mm (+4 mm) are depicted, plus the tibial metal baseplate fixed to the resected proximal tibia.

To allow for an investigation of the model kinetic predictions on the TF joint, the primary outcome measures of the simulations included the medial collateral ligament (MCL), lateral collateral ligament (LCL), and posterior cruciate ligament (PCL) forces and strain, and the TF compressive force. To assess the biomechanical effects on the PF joint, the secondary simulation outcomes were the forces around the knee extensor mechanism, including the quadriceps muscle force, quadriceps tendon-to-femur force, and PF compressive force. The lateral PF ligament bundles were found to remain slack, while the medial PF ligament forces were essentially identical throughout the series of trial simulations regardless of the PE insert thickness and hence excluded from the subsequent analysis. In addition, the PF joint kinematics, including patellar flexion, anterior translation, distal translation, medial rotation, lateral tilt, and lateral shift, were recorded. To quantify the biomechanical effects of PE insert thickness across the joints, we identified the average and peak values in the ligament forces, TF and PF compressive forces throughout the simulated range of motion; and subsequently computed the average and peak percentual change with respect to the reference case.

#### **3. Results**

#### *3.1. Tibiofemoral Joint*

The peak ligament and compressive forces on the TF joint for all the simulated PE insert thickness cases are depicted in Figure 3. The peak MCL, LCL, and PCL forces increased by 38 N, 74 N, and 125 N, respectively, as PE insert thickness changed from 9 mm to 11 mm; the corresponding peak force increase for a 4 mm greater thickness relative to the reference case was 80 N, 157 N, and 286 N. Conversely, the peak force in the MCL, LCL, and PCL decreased by 36 N, 52 N, and 98 N, respectively, when using a 2 mm thinner PE insert compared to the reference case, while with a 4 mm smaller thickness than the reference, the peak decrease was 57 N, 70 N, 177 N, respectively. Table 1 summarizes the peak as well the average percentage variations relative to the reference case. Varying PE insert thickness simultaneously changed the ligament individual bundle strain. Both medial and lateral collateral ligaments appeared to stretch in the same strain range of about 0–6% for all the simulated thickness cases, whereas the PCL exhibited a maximum strain around 10% at 90 ◦ of knee flexion with 13 mm thickness. Force and strain patterns of the MCL, LCL, and PCL at varying thicknesses are provided in Appendix A.

**Figure 3.** Peak ligament and compressive forces on the tibiofemoral (TF) joint at varying PE insert thicknesses during a chair-rising simulation. From left to right: MCL force, LCL force, PCL force, and TF compressive force. Results are displayed relative to the reference case (gray filled marker).

**Table 1.** Changes in muscle, ligament, and joint compressive forces due to tibial PE insert thickness variations.


Data are presented as average and peak percentage difference over the chair-rising simulation for each 2 mm PE insert thickness change relative to the reference case.

> Changes in PE insert thickness affected the TF compressive force, denoting a larger effect at 90 ◦ of knee flexion angle at which the highest peak occurred at 4.5 times body weight (BW) with a 13 mm PE insert (Figure 3). Compared to the reference case, the peak TF compressive force increased by 0.3 BW and 0.8 BW with 2 mm and 4 mm greater PE insert thicknesses, respectively. Using 2 mm and 4 mm thinner PE inserts than in the reference case showed a peak decrease of 0.2 BW and 0.3 BW, respectively. Detailed peak and average percentage variations of the TF joint compressive force are also provided in Table 1.

#### *3.2. Patellofemoral Joint*

The forces around the knee extensor mechanism exhibited peaks during the momentum transfer phase of the chair-rising movement at about 90 ◦ of knee flexion. Peak values of the quadriceps muscle force, quadriceps tendon-to-femur force, and PF compressive force for all the simulated PE insert thickness cases are shown in Figure 4, and their variations with regard to the reference case are summarized in Table 1. Interestingly, the PF compressive force followed an opposite trend than the compressive force on the TF side, indicating a peak decrease of 0.1 BW and 0.3 BW when PE insert thickness was varied from 9 mm to 11 mm or 13 mm, respectively. In contrast, simulation of thinner tibial PE inserts, relative to the reference case, resulted in a slightly increased peak PF joint compressive force by 0.1 BW (−2 mm) and 0.3 BW (−4 mm). The force–angle curves of the knee extensor mechanism structures are provided in Appendix A.

**Figure 4.** Peak forces around the knee extensor mechanism at varying PE insert thicknesses during a chair-rising simulation. From left to right: quadriceps muscle force, quadriceps tendon-to-femur force, and patellofemoral (PF) compressive force.

Changes in PE insert thickness affected the patellar proximal–distal translation and medial–lateral rotation, with a minimal effect at full knee flexion (Figure 5). The patella shifted by 0.8 mm distally for every 2 mm incremental change of the PE insert thickness, from 9 mm (reference) up to 13 mm, and about 0.8 more proximally for each 2 mm decrement relative to the reference thickness. Compared to the reference case, the patella was rotated 0.4 and 0.7 mm more medially with 11 mm and 13 mm PE insert thicknesses, respectively and, accordingly, 0.4 and 0.1 mm more laterally with a 7 mm and 5 mm PE insert thickness, respectively.

**Figure 5.** Reference frames used to express the PF kinematics (left). Muscles and ligaments of the TF and PF joints are hidden in this model view. Patellar kinematic profiles at varying PE insert thicknesses during a chair-rising simulation (right).

#### **4. Discussion**

The most important findings in the present study were that a larger thickness of the tibial PE insert elevated the ligament forces throughout the range of knee flexion and extension and, consequently, increased the TF compressive force across the mediolateral compartment; in contrast, varying the PE insert thickness had little effect on the PF biomechanics, although it indicated an unexpected slight decrease in the PF joint loading. These results confirm our hypothesis and suggest that there is a trade-off in the kinetic behavior between the TF and PF joint structures with regard to varying PE insert thickness.

Tibial PE insert thickness variation had a marked effect on the ligament forces and strain patterns, as expected because of overstuffing the TF joint space. With thicker PE inserts relative to the reference case, the MCL, LCL, and PCL forces increased considerably both in flexion and extension, although the changes in the PCL were more distinct after mid-flexion. The observed higher LCL forces over the MCL might be related to the larger stiffness value assigned to the MCL individual bundles (3000 N) based on literature experimental evidence [26]. Another fact worth noticing is that the PCL was engaged at lower flexion angles with thicker PE inserts, which may be due to the larger joint distraction in extension. On the other hand, thinner PE inserts (5 mm, 7 mm) slackened the MCL and LCL almost entirely from 60◦ to 90◦ , due to the larger flexion and extension gap. This is an unfavorable scenario, as slackening the collateral ligaments may destabilize the knee in flexion. The PCL, unlike the other two ligaments, received more tension in the flexed knee with PE inserts thicker than 9 mm, reaching a peak strain at about 10% when using a 13 mm component. This is, however, well below the yield strain of 14% as reported in the literature, above which there might be a ligament injury as a result of overstretching [22]. Similar to the ligaments' kinetic behavior, the compressive force in the TF joint increased substantially at varying PE insert thicknesses, particularly at 90◦ of knee flexion with a 13 mm PE component. Considering that the quadriceps muscle force was roughly the same across PE insert thickness cases, this increase can be solely attributed to the summation of increased joint ligament forces. The magnitude of the TF compressive force found in this study is consistent with earlier observations [27], reporting an average peak of about four times BW during rising from a chair without the aid of arms in a natural knee. This supports our reported values with 9 mm (3.7 BW) or 11 mm (4 BW) thicknesses and raises the concern that using thicker PE inserts, such as 13 mm (4.5 BW) or more, could greatly elevate the joint loads on the PE surface, leading to destructive wear [28]. Surgeons should consider this surgical option relative to the resected bone on the tibial plateau, bearing in mind that a thicker PE insert also requires a lower level of tibial resection, which is associated with posteromedial bone failure and early aseptic loosening [29]. In addition, deeper tibial resections may lead to the removal of a considerable part of the tibial PCL attachment, which can result in a reduced femoral rollback and anteroposterior instability during flexion [30]. A recent cadaveric study proposed a PE insert thickness of 10 mm with a posterior slope of 5◦ to preserve more than 50% of the tibial PCL attachment site [31]. Our findings highlight the importance of preoperative planning of the appropriate tibial insert thickness, which considers both the thickness value available and the tibial resection depth and, additionally, support the conceptual premise that selection of a tibial insert that is either too thin or too thick can have a negative impact on the TF kinetics, leading to a sub-optimal solution.

Varying the thickness of the PE insert revealed marginal secondary effects on the PF biomechanics. The quadriceps muscle force remained unchanged across PE insert thickness variations in the full arc of the knee range of motion. However, more thickness in the tibial PE component increased the quadriceps tendon-to-femur force in the range of 60◦–90◦ of knee flexion. Contrary to our expectations, the PF compressive force slightly decreased with 2 mm and 4 mm greater PE insert thickness, compared to the reference case, and this change was more perceivable from mid-flexion to 90◦ . This is in agreement with the findings of a previous study [32] and might be explained by quadriceps–femur load sharing as the quadriceps wraps around the distal femur in a flexed knee position. This mechanism is also reflected in the PF kinematics. The patella was consistently more distal with thicker PE inserts than in the reference case throughout the range of flexion and extension, as expected due to the joint line elevation, which subsequently results in patella baja. From a clinical point of view, unloading the PF joint could reduce patellar complications after TKA, especially anterior knee pain, which strongly correlates with postoperative patient

dissatisfaction and impaired quality of life [33–35]. On the other hand, patella baja may have serious consequences on the overall function of the knee after surgery, such as limited range of motion because of patellar maltracking. Moreover, the distally displaced patella can impinge on the anterior part of the tibial PE insert or the tibial tray during flexion, potentially increasing wear [36]. Hence, surgeons should carefully assess the concomitant effects on the PF joint when pre-planning the desired PE insert thickness to ensure that the patella slides properly in the trochlear groove.

To the best of our knowledge, this is the first musculoskeletal simulation study which explores the biomechanical consequences of tibial PE insert thickness variations on both the TF and PF joints after cruciate-retaining TKA. A major strength of this study is the use of a computational modeling approach to investigate the effect of tibial PE insert thickness in a highly controlled simulation environment, where all the other surgical variables, such as the tibial resection depth, remained unchanged. This overcomes an important limitation of clinical studies, in which confounding factors, such as the implant design, size, and alignment, are present and can potentially affect the surgical outcomes [9].

This study has several limitations. At first, the modeled ligament mechanical properties, such as stiffness and reference strain, were determined from the literature. Furthermore, the musculoskeletal simulation represented only a mechanically aligned cruciateretaining TKA, meaning that the results cannot be extrapolated to other surgical techniques or implant designs. It is also noteworthy that the computational model simulated only one patient, disregarding the anatomical variability or other pathological conditions among different patients. Further research to overcome this limitation may be to characterize the effect of variability in the morphological knee joint phenotypes on the postoperative kinetics with varying tibial PE insert thicknesses. A promising methodology to explore this could be the combination of statistical shape modeling and musculoskeletal simulation.

#### **5. Conclusions**

Changes in tibial PE insert thickness have considerable influence on the kinetics across the TF articulation and surrounding ligamentous structures, but marginal effects on the PF biomechanics. Increasing the tibial PE insert thickness resulted in higher ligament forces, and subsequently increased loading across the medial and lateral TF compartments in the full arc of knee motion after cruciate-retaining TKA. However, thicker PE inserts resulted in a slightly lower PF force by means of an increased load-sharing between the quadriceps tendon and femur.

**Author Contributions:** Conceptualization, P.T. and N.V.; methodology, P.T. and M.A.M.; software, P.T., M.A.M., and R.F.; validation, P.T., M.A.M., and R.F.; formal analysis, P.T.; investigation, P.T.; resources, P.T. and M.A.M.; data curation, P.T., M.A.M., and R.F.; writing—original draft preparation, P.T.; writing—review and editing, P.T., M.A.M., R.F., and N.V.; visualization, P.T.; supervision, M.A.M., N.V., and B.K.; project administration, N.V. and B.K.; funding acquisition, M.A.M. and N.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Top Technology Twente Connecting Industry program (TKI Top Sector HTSM) and Stryker European Operations Ltd., Ireland.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We would like to thank Eric Garling (Stryker, Kalamazoo, MI, USA) and Arman Motesharei (Stryker, Newbury, UK) for providing the geometry files of Triathlon prosthesis used in this study.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Figure A1.** Ligament force and strain patterns at varying PE insert thicknesses during a chair-rising simulation. From left to right top and bottom: MCL force and strain, LCL force and strain, PCL force and strain.

**Figure A2.** TF compressive force at varying PE insert thicknesses during a chair-rising simulation.

**Figure A3.** Forces around the knee extensor mechanism at varying PE insert thicknesses during a chair-rising simulation. From left to right: quadriceps muscle force, quadriceps tendon-to-femur force, and PF compressive force.

#### **References**


## *Article* **Effects of Hip Abductor Strengthening on Musculoskeletal Loading in Hip Dysplasia Patients after Total Hip Replacement**

**Giordano Valente 1,\* , Fulvia Taddei <sup>1</sup> , Alberto Leardini <sup>2</sup> and Maria Grazia Benedetti <sup>3</sup>**


**Abstract:** Hip dysplasia patients after total hip replacement show worse functional performance compared to primary osteoarthritis patients, and unfortunately there is no research on muscle and joint loads that would help understand rehabilitation effects, motor dysfunctions and failure events. We tested the hypothesis that a higher functional improvement in hip dysplasia patients who received hip abductor strengthening after hip replacement, would result in different gait function and musculoskeletal loads during walking compared to patients who performed standard rehabilitation only. In vivo gait analysis and musculoskeletal modeling were used to analyze the differences in gait parameters and hip and muscle forces during walking between the two groups of patients. We found that, in a functional scenario of very mild abnormalities, the patients who performed muscle strengthening expressed a more physiological force pattern and a generally greater force in the operated limb, although statistically significant in limited portions of the gait cycle, and likely related to a higher gait speed. We conclude that in a low-demand task, the abductor strengthening program does not have a marked effect on hip loads, and further studies on hip dysplasia patients would help clarify the effect of muscle strengthening on loads.

**Keywords:** total hip replacement; hip dysplasia; hip abductor muscles; musculoskeletal modeling; gait analysis; rehabilitation

#### **1. Introduction**

Developmental dysplasia of the hip (DDH) features chronic dislocation of the femoral head, leading to secondary osteoarthritis, contractures and atrophy of the hip abductor muscles, and resulting in overall motor dysfunction, pain and disability [1]. These musculoskeletal issues can be treated through total hip replacement (THR), which includes a specific complexity for DDH patients due to difficulties in restoring the ideal center of joint rotation. The major risks after THR are residual abnormalities in gait and early wear of the prosthesis, due to poor abduction function and femoral offset [2,3]. Rehabilitation of the hip abductor muscles has been debated and targeted as a major player in the removal of residual dysfunctions [4,5], although no consensus has been found on a specific rehabilitation protocol for hip abductor strengthening [6,7]. Recent studies on THR patients following osteoarthritis demonstrated the effectiveness of hip muscle strengthening and activity retraining intervention on functional outcomes [8,9].

Recently, we have evaluated the effects of a program of abductor exercises on a population of DDH patients through a randomized prospective study [10]. We found a significantly higher improvement in clinical and functional scores at follow-ups of a Study Group that performed a specific muscle strengthening program, compared to a Control Group that performed standard postoperative rehabilitation only. We concluded that a

**Citation:** Valente, G.; Taddei, F.; Leardini, A.; Benedetti, M.G. Effects of Hip Abductor Strengthening on Musculoskeletal Loading in Hip Dysplasia Patients after Total Hip Replacement. *Appl. Sci.* **2021**, *11*, 2123. https://doi.org/10.3390/ app11052123

Academic Editor: Carlo Albino Frigo

Received: 18 December 2020 Accepted: 24 February 2021 Published: 27 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

specific protocol for hip abductor strengthening with adequate presurgical planning would improve functional performance and patient satisfaction in DDH patients after THR.

From a biomechanical perspective, gait analysis studies among THR patients demonstrated that DDH patients have worse functional performance and unsatisfactory clinical results compared to primary osteoarthritis patients [11,12]. However, although a few studies analyzed the effect of hip abductor strengthening on joint loads in osteoarthritis patients after THR [13,14], there is no literature regarding muscle and joint contact forces during movement in DDH patients after THR. Muscle and joint loads are related to the stress at the hip joint, with consequent failure events that may occur as a result of wear, edge loading, impingement and subluxation [15,16].

Quantifying muscle and joint loads after rehabilitation following THR of DDH patients would provide insights into the biomechanical effect of a program of exercises, identify potential compensation strategies, or help explain pain and residual motor dysfunctions. The state-of-the-art method to quantify musculoskeletal loads includes the use of multibody models of the musculoskeletal system and patient-specific gait analysis measurements, with consequent simulations of movement and prediction of the biomechanical variables [17,18].

Therefore, the aim of this study was to analyze the effect of a rehabilitation protocol for abductor muscle strengthening on gait function and musculoskeletal loads during walking in DDH patients after THR. We compared gait analysis measurements, joint contact forces and muscle forces of the hip between a group of patients who performed the specific rehabilitation protocol and another group of patients who performed standard rehabilitation only, in the operated and contralateral limb. We tested the hypothesis that a higher improvement in clinical and functional outcome in patients after hip abductor strengthening, as shown recently [10], would result in differences in muscle and joint forces compared to patients who performed standard rehabilitation only.

#### **2. Materials and Methods**

#### *2.1. Patients and Experimental Data*

Thirty DDH patients were randomly selected from a total of 110 recruited previously within a multicentric clinical trial [10] about the effects of a strengthening abductor program in functional performance and patient's satisfaction. The 30 patients were scheduled for gait analysis at a follow-up of at least 6 months after THR [19–21]. Ten of them did not complete the protocol refusing to perform gait analysis, and of the 20 patients that accepted to perform the exam, 13 were in the rehabilitation arm (Study Group) and seven in the control arm (Control Group) of the study. After gait analysis data acquisition, one patient in the Study Group and one patient in the Control Group were excluded because of missing data on ground reaction forces, therefore 18 patients were included.

All patients were operated with a minimally invasive anterolateral approach for osteoarthritis secondary to DDH [10], and treated with the same cementless hip implant. Preoperative planning was carried out with HipOp-Plan software (http://www.hipopplan.com/), which allows patient-specific, interactive and user-friendly positioning of the prosthesis components.

Both groups underwent standard early postoperative rehabilitation based on mobilization, gait training and generic muscle strengthening, until patient discharge. In addition, the Study Group only, performed a 2-week daily rehabilitation program aimed at strengthening the hip abductor muscles [8]. Measurements such as Western Ontario and McMaster Universities Osteoarthritis Index (WOMAC), Harris Hip Score and 10 Meter Walk Test were used to measure the functional outcome.

The institutional Ethics Committee approved the study (n. 0006168, registered in the local Institutional Registry with the number 26/13), and all patients read and signed an informed consent.

#### *2.2. Gait Analysis*

After completion of the rehabilitation program, the patients performed state-of-the-art gait analysis. Each subject performed five trials of barefoot level walking at a self-selected speed on a 10 m walkway. An additional static trial was recorded in standing upright position. The subjects were first instrumented with 22 reflective markers on pelvis and lower limbs following an established protocol [22]. Three-dimensional marker trajectories and ground reaction force data were collected using an 8-camera motion capture system (100 Hz, Vicon 612 Motion System, Oxford, UK), and two embedded force platforms (2000 Hz, Kistler, Winterthur, Switzerland), respectively.

#### *2.3. Musculoskeletal Modeling and Gait Simulations*

A generic multibody musculoskeletal model including 14 rigid bodies, 20 degrees of freedom and 80 musculotendon actuators [23] was used in conjunction with experimental marker trajectories and ground reaction forces to calculate joint angles, joint moments, muscle forces and joint contact forces during walking, by leveraging OpenSim [24]. Preliminarily, marker trajectories were processed with a Woltring low-pass filter to smooth the data with a mean square error of 10 mm, then MATLAB (R2018a, Mathworks, USA) was used to process marker trajectories and ground reaction forces to generate the input files and batch process the simulations in OpenSim. First, the model was personalized for each subject, by scaling the dimensions of each body segment, mass and inertial properties, and the elements attached to the body segments, based on (i) the distances between the experimental markers from the static trial and the corresponding virtual markers on the model, and (ii) body mass. Joint angles during walking were then calculated through inverse kinematics, by minimizing the errors between experimental and virtual markers. Joint moments were calculated through inverse dynamics, from the joint angles and the ground reaction forces. Then, muscle forces were calculated by decomposing the joint moments among the musculotendon actuators through static optimization, by minimizing the sum of muscle activation squared and accounting for the force–length–velocity relationship [25]. Finally, joint contact forces were calculated from the instantaneous force equilibrium through joint reaction analysis.

All data were normalized to percentage of gait cycle from 1% to 100%, and muscle and joint contact forces were expressed in body weight (BW). Hip contact forces, and the most influential muscle forces acting on the hip were analyzed. The muscle forces were grouped and summed in hip abductor muscles (gluteus medius anterior, middle and posterior, gluteus minimum anterior, middle and posterior) and hip flexor muscles (iliacus, psoas and rectus femoris).

#### *2.4. Statistical Analyses*

Differences in subject characteristics, functional parameters and peak values in gait analysis measurements between the Study Group and the Control Group were tested by using a Mann–Whitney U-test (α = 0.05) for independent nonparametric samples.

The differences in hip contact forces and muscle forces between groups across the gait cycle were analyzed in both the operated and the contralateral limb. First, the forces were expressed as mean and standard deviation, and mean differences between groups were calculated to have a preliminary quantification of the differences for an initial overall comparison. Then, to evaluate group-wise statistically significant differences, statistical parametric mapping (SPM) was used [26]. Specifically, nonparametric two-tailed unpaired t-tests were performed between the variables of the Study Group and the Control Group (α = 0.05), by using the SPM1D package (SPM, www.spm1d.org, v0.4 [27]) in MATLAB. Differences were considered clinically relevant if significant differences occurred for at least 4% of the gait cycle [28].

Further, an analysis on the overall peaks of hip contact forces was performed. First, force peaks across the gait cycle were calculated and presented as boxplot distributions, and statistically significant differences between groups were tested by using a Mann–Whitney

U-test (α = 0.05). Finally, linear correlation (R, *p*-values) between force peaks and gait speed was checked in both the Study Group and the Control Group patients.

#### **3. Results**

#### *3.1. Function and Gait Analysis Results*

The general characteristics of the patients included are reported in Table 1. Age, body mass, height, and BMI showed no significant differences between the two groups.

**Table 1.** General characteristics of the patients (Median (IQR)).


Measurements from gait analysis demonstrated that cadence and gait speed were significantly higher in the Study Group, while stance time in the operated side and cycle time were significantly higher in the Control Group (Table 2). The other time–distance parameters showed no significant differences.

**Table 2.** Comparison of time–distance parameters.


Regarding kinematic measurements, there were no significant differences in pelvic obliquity and hip range of motion between the Study Group and the Control Group (Table 3).


**Table 3.** Comparison of maximum peaks in pelvic and hip kinematics.

#### *3.2. Analysis of the Musculoskeletal Forces*

Regarding the forces predicted via musculoskeletal modeling, the hip contact forces in the operated leg were larger in the Study Group throughout most of the gait cycle, with a maximum difference of 0.6 BW, although the differences were significant in limited portions of the gait cycle (Figure 1).

Conversely, hip contact forces in the contralateral leg were larger in the Control Group, with a maximum difference of −0.6 BW, although not statistically significant (Figure 1).

However, from a qualitative analysis in the shape of the hip contact forces, the Study Group shows a more physiological pattern where a typical double-bump is present with a central absorption loading phase, in both the operated and contralateral leg. No remarkable differences were found in the hip abductor and flexor muscles.

Focusing on the peaks of hip contact forces, no significant differences between the Study and the Control Group were found in the operated limb, while significant differences were present in the contralateral limb with higher force peaks in the Control Group (Figure 2).

In addition, there was a moderate significant correlation between gait speed and contact force peaks in the hip in the Study Group (R = 0.44, *p* < 0.001) (Figure 3). The peak of hip contact force occurred always toward toe off.

**Figure 1.** Hip contact forces and muscle forces (normalized to body weight and represented as mean and standard deviation) of the operated leg (**A**) and contralateral leg (**B**) across the gait cycle, compared between the Study Group and the Control Group. The analysis is focused on the hip, including joint contact forces and major muscle forces grouped by function. The asterisk (\*) indicates the maximum mean difference. Results from SPM t-tests are reported: when SPM t-values exceed the critical threshold (horizontal dotted lines), the differences between groups are significant, and when significant differences occurred for at least 4% of the gait cycle, the areas are reported in grey with the corresponding *p*-values.

**Figure 2.** Boxplot distributions of the peaks of hip contact forces during gait compared between the Study and Control Groups, in the operated leg (**A**) and contralateral leg (**B**). *P*-values are reported for the significant differences between groups (Mann–Whitney U-test).

**Figure 3.** Linear correlations between gait speed and hip contact force peaks in all walking trials, for the Study Group (**A**) and Control Group (**B**). Correlation coefficients (R) are reported and indicated with an asterisk (\*) when statistically significant (*p* < 0.05).

#### **4. Discussion**

In this study, we used in vivo motion analysis and multibody musculoskeletal models to test the hypothesis that a higher improvement in clinical and functional outcome shown in DDH patients who received hip abductor strengthening after hip replacement [10], would result in different musculoskeletal loads (i.e., joint contact forces and muscle forces) during walking compared to patients who performed standard rehabilitation only.

Gait analysis revealed mild abnormalities in the time–distance parameters and in the pelvic and hip kinematics, with few significant differences between the two groups (Tables 2 and 3). Although only few articles are available for gait analysis in DDH patients operated on hip prosthesis [11,19,29], our results are in general agreement with findings of slower gait, shorter stride, reduced range of motion in the sagittal and coronal planes [11], and lower pelvic rotation in the coronal plane [29] with respect to healthy subjects. In

addition, functional outcome data in terms of 10 Meter Walk Test, Harris Hip Score and WOMAC demonstrated a better outcome in the Study Group with results very similar to those of the entire sample of patients reported in [10] (Table S1).

Overall, the protocol for abductor muscle strengthening has a mild effect on joint contact forces and muscle forces during walking in the operated limb. We would have expected a more marked effect of strengthening on loads from the hypothesis tested, given the significant difference in abductor strength between groups [10], and the major influence of hip abductor strength on joint contact forces shown in a previous simulation study [30]. In the present study, the patients in the Study Group generally experienced a greater hip force in the operated limb compared to the Control Group (Figures 1 and 2), although with significant differences in a small portion of the gait cycle. This mild difference in hip contact force is, however, likely related to gait speed, which was greater in the Study Group (Table 3, Figure 3). The influence of gait speed on hip contact force during the gait cycle can also explain the force pattern where a typical two-peak pattern is present, more similar to that in both healthy population and THR patients, with a lower force magnitude in both peaks compared to the healthy [28,31]. As we found speed as an influent parameter in our results (Figure 3), we performed a post-hoc analysis where we tried to remove the two slower walking patients from the six patients of the Control Group, with all possible cautions in interpreting the results, given the small number of subjects. This led to a more physiological pattern of hip contact force (i.e., more marked double-bump), but the statistical significance of the differences between groups was only changed slightly. Therefore, we decided to keep all patients in the Control Group, considering also that the two slower patients were globally not more nor less involved than the others.

The lack of differences in muscle forces of the hip abductors and flexors between limbs is difficult to interpret, also considering that there is no consensus in the literature about the effect of impaired muscle status on joint contact forces.

Unfortunately, no similar studies are available on hip contact and muscle forces after THR in DDH patients, indeed this is the first study that analyzed the musculoskeletal loads during gait in DDH patients after THR. Previous research focused on dysplastic hips [32,33] or in patients operated on periacetabular osteotomy [34,35]. Sørensen et al [35] reported a normal hip muscle function during walking after 12 month periacetabular osteotomy, and found that joint force magnitude continued to be higher than normal, justifying this finding with the need for joint structures of a longer time to heal than muscles, with a residual pain alleviating strategy. Damm et al. [36] reported a correlation between muscle atrophy, fatty degeneration and higher hip joint forces. Fatty atrophy was observed after THR contributing to functional outcome deterioration [37].

Considering that the patients included in the present study underwent CT-based planning that allowed targeting more accurately the reconstruction of the center of rotation and other features of the implant, we can suppose that the hip contact force in the Study Group could represent an optimal outcome in terms of recovery of correct loading, also symmetric with the contralateral hip. Nevertheless, in the Control Group, the hip force peak seems to be lower, although not significantly, in the operated limb and significantly higher in the contralateral limb, which suggests a nonoptimal recovery of the operated hip function. According to Sørensen et al. [35], who found contradictory results assessing hip joint kinetics in DDH patients by isokinetic dynamometry and gait analysis, walking can be considered a low-demand task in terms of muscle recruitment, thus more strenuous activities should be explored in the future for a more robust functional assessment of hip muscles. Actually, the physical demand on the gluteal muscles is higher during other common motor tasks (e.g., stairs, squat), and the effect of a more performing muscle might be more marked [14,36,38,39].

The main limitation of this study is the small sample size of the patients, with relevant differences between the Study Group and the Control Group. Unfortunately, we were not able to control the adherence of the patients to the original protocol after 6 months from surgery. In fact, 10 of the patients randomly included in the subgroup for gait analysis came from a distance, and at 6-month follow up they declined to participate. This limited sample size has undoubtedly limited the significance of the observed differences. However, we believe that this limitation does not invalidate the study, given the absence in the literature of analogous studies and the importance of the topic, as it could bring preliminary results to the discussion. Another limitation might regard the lack of gait analysis data before surgery. It is indeed true that a patient control pre- and postoperation would be of additional help in interpreting the results, but the gait pattern before surgery is usually confounding for the presence of leg length discrepancy and pain, hence this choice might be questionable anyway. Finally, it became evident from the analysis of the results that the chosen motor task (i.e., walking at self-selected speed), besides being important as the most common daily task, might not be the best task to highlight the differences induced in kinematics, joint and muscle loads between the operated and contralateral limb. It is possible that a higher demanding task, such as squat or stairs, might better highlight the differences. This should be kept in mind when designing similar studies in the future.

#### **5. Conclusions**

In conclusion, in this pilot study we found that the protocol for abductor muscle strengthening in hip dysplasia patients following total hip replacement does not have a major effect on muscle and joint contact forces during walking. However, in a functional scenario of very mild abnormalities, patients who performed muscle strengthening expressed a more physiological force pattern and a generally greater force in the operated limb, both likely related to a higher gait speed. Therefore, in a low-demand task such as walking, the abductor strengthening program leads to more physiological load patterns, although the differences in magnitude are not marked. Further studies on a larger sample size of DDH patients are desirable to clarify the effects of abductor muscle strength on hip joint loads.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2076-341 7/11/5/2123/s1, Table S1: Functional outcome measurements.

**Author Contributions:** Conceptualization: G.V., F.T., M.G.B.; data acquisition: A.L.; modeling and simulations: G.V.; data analysis: G.V., M.G.B.; writing—original draft preparation: G.V.; writing review and editing: G.V., F.T., A.L., M.G.B.; funding acquisition: M.G.B. All authors have read and agreed to the present version of the manuscript.

**Funding:** This study was funded by the Programma di ricerca Regione-Università 2010–2012 Area 2—"Ricerca per il Governo clinico", within the project title "Optimising functional recovery in patients with developmental dysplasia of the hip: an innovative pathway for joint reconstruction and rehabilitation".

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of the Rizzoli Orthopaedic Institute (protocol n. 0006168, registered in the local Institutional Registry with the number 26/13).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


forces in patients with hip dysplasia. *Gait Posture* **2015**, *42*, 529–533. [CrossRef]


*Article*

## **The E**ff**ects of the Rectus Femoris Muscle on Knee and Foot Kinematics during the Swing Phase of Normal Walking**

#### **Carlo Albino Frigo 1,\* , Christian Wyss 2,3 and Reinald Brunner 2,3**


Received: 10 September 2020; Accepted: 5 November 2020; Published: 6 November 2020 -

**Abstract:** The role of rectus femoris (RF) muscle during walking was analyzed through musculoskeletal models to understand the effects of muscle weakness and hyperactivity. Such understanding is fundamental when dealing with pathological gait, but the contribution of RF as a bi-articular muscle is particularly difficult to estimate. Anybody software was used for inverse dynamics computation, and SimWise-4D for forward dynamics simulations. RF force was changed in the range of 0 to 150%, and the resulting kinematics were analyzed. Inverse dynamics showed a short positive RF power in correspondence with the onset of knee extension in the swing phase. Forward dynamics simulations showed an increasing knee flexion and initial toe contact when the RF force was decreased, and increasing knee extension and difficult foot clearance when the RF force was increased. The step became shorter with both increased and reduced RF force. In conclusion, the RF actively contributes to the knee extension in the swing phase. RF also contributes to obtaining a proper step length and to preparing the foot for initial heel contact. So the effect of RF muscle as a bi-articular muscle seems fundamental in controlling the motion of distal segments. RF overactivity should be considered as a possible cause for poor foot clearance in some clinical cases, while RF weakness should be considered in cases with apparent equinus.

**Keywords:** rectus femoris; foot kinematics; knee kinematics; normal gait; gait cycle; forward modeling; musculoskeletal modeling

#### **1. Introduction**

Although the function of a muscle is determined by its anatomical arrangement, the real effect of muscle contraction during function depends on the state of contraction of other muscles, on the external forces applied to the body, and on the so-called "dynamic coupling" among the anatomical segments. For example, during walking, excessive activity of the ankle plantar flexors at the beginning of the stance phase can lead to knee hyperextension and pelvis internal rotation in cerebral palsy children [1]. A similar contraction occurring at the late stance phase can produce external rotation of the pelvis and, depending on the level of gastrocnemius contribution, knee flexion instead of knee hyperextension. Quadriceps and hamstring muscle co-contractions can have different effects at different joints configurations [2], and it appeared that the activity of vasti (VA) associated with hamstrings can enhance the extensor effect of these muscles at the hip, particularly when the knee and the hip are relatively flexed.

To understand the role of the muscles involved in locomotion a study should include an analysis of the muscle activity and also an estimate of what happens if the muscle activity changes.

In the present work, we adopted this approach to investigate the role of rectus femoris (RF) on the execution of the swing phase during walking. In fact, several factors acting at late stance-early swing can affect the swing phase: hip acceleration, ankle plantar flexion, ground reaction forces, inertia, but the ability of RF (a proximal, double joint muscle) to control the distal segments seems to play a relevant role in controlling the movement amplitude and the position of knee and foot at the end of the swing phase. From a clinical point of view, a more thorough understanding of the effects of the RF activity can be useful in planning treatment. Swing phase pathologies are difficult to understand, especially because the knowledge on the effect of RF during swing already in normal conditions is insufficient. If it changes activity, the mechanical effects can be even less well estimated.

Previous studies have shown that the RF activity is present at late stance-early swing to control knee flexion [3–5]. Nene et al. [6] showed that the pattern of this activity depends on gait velocity and also on the angular velocity of the shank. Although some crosstalk of electromyographic (EMG) signals occurs between VA and RF [7], the same authors have concluded that during the initial swing the recorded activity mostly originated from the RF. Additionally, in a recent publication [8], the most recurrent pattern of activation was shown to be composed of three phasic activities: at the beginning of the gait cycle, around toe-off, and at the terminal swing. The activity corresponding to the late stance-early swing was considered consistent with the need to control the knee flexion. Muscle-actuated forward dynamic simulations of the swing phase of normal gait further demonstrated that the RF is important to control knee flexion during swing [9]. Its removal resulted in increased knee flexion during the swing phase which was even more pronounced in pre-swing than in early swing [10]. According to other authors [11], the RF was mainly active when a combination of both hip flexion and knee extension moments were required, and not when the knee extension moment was combined with a hip extension moment. This is consistent with the anatomical arrangement of the RF that acts as a hip flexor and knee extensor, and thus it seems particularly suitable when both of these effects are required during movement. However, the following questions remain open and this study was aimed at finding a fundamental answer to them:

Does the RF only control knee flexion in the early swing phase by breaking the movement or does it also actively extend the knee in the second phase?

To what extent are joint kinematics and the segmental position at initial contact altered when the RF is removed?

To what extent are joint kinematics influenced by RF hyper-activity at late stance-early swing?

#### **2. Materials and Methods**

RF activity during gait has been described by several authors [3–8]. We in contrast, were not interested in activity patterns, but in data of one single normal individual which we modified in a forward modelling study in order to understand the principle of the RF effect. Difference in segment sizes and proportions depending on age may change the values but not the principle as the gait pattern after four years of age changes little [12].

(a) Data acquisition and models.

We performed an inverse dynamics analysis as a first step. Since our study was mostly theoretical, we just used an exemplary collection of data from one normal subject as the input to our model. Kinematic and kinetic data referred to a 23-years-old woman (weight 59.1 kg; height 164 cm) who was walking barefoot at her natural speed in a gait analysis laboratory (motion analyzer Vicon Nexus 1.8.5, 12 cameras, Oxford Metrics plc, Oxford, United Kingdom; 4 dynamometric platforms, Kistler, Winthertur, Switzerland; and a 12 channels Noraxon EMG system, Noraxon, Scottsdale, AZ, USA). Data were processed according to a standard protocol [13], and were imported into a musculoskeletal model of the AnyBody Software Version 7.0 (AnyBody Technologies, Aalborg, Denmark). An average

of six trials was used for the inverse dynamic simulation, which allowed us to obtain joint moments and individual muscle forces and muscle power.

In a second phase, normal gait data were imported into a software platform (SimWise-4D, Design Simulation Technologies, Canton, MI, USA) in which a model of a human subject was constructed for dynamics simulation, with some additional features with respect to a previously developed model [1,2]. The model was composed of solid segments representing the trunk, the pelvis, the thigh, the shank, and the foot for both the left and the right sides. The mass and moments of inertia of these solids were tailored on the subject based on anthropometric tables [14]. The segments were connected to each other by revolute joints, allowing several degrees of freedom: three between trunk and pelvis, three at the hips, two at the knees (flexion/extension, internal/external rotation), and two at the ankle (dorsi/plantar flexion, inversion/eversion). Each revolute joint could be controlled by the corresponding joint angle time course, so that the normal kinematics could be reproduced. In this case the joint moments resulting from the dynamic equilibrium of forces and moments applied to the body segments were obtained (inverse dynamics problem). In a different approach, the revolute joints could be left free to rotate (no internal moment was produced) and the model provided the kinematics resulting from forces and moments applied to the system (forward dynamics problem). Digital anatomical models of pelvis, femur, tibia, and fibula were superimposed to the corresponding solid segments to obtain a reference for the muscle attachment points. Several muscles of the lower limb, represented by linear actuators, were then attached to the model: gluteus maximus, iliopsoas, semimembranosus, semitendinosus, biceps femoris long head and short head, gastrocnemius medialis and lateralis, soleus, tibialis anterior, RF, and VA (intermedius, lateral, medial). The distal tendons of some muscles (gluteus maximus, iliopsoas, semimembranosus, quadriceps femoris, gastrocnemii) were reproduced by a chain of small cylinders that could slide over wrap surfaces representing the bone surfaces. The quadriceps tendon was connected to a parallelepiped representing the patella. This, in turn, was connected to the anterior tuberosity of the tibia by a rope representing the patellar tendon. The patella could slide over a wrap surface representing the femoral trochlea. For the present use of the model, all muscles except the RF were considered inactive (just a 1 N force was imposed to keep the tendons in place). The RF force instead, was input into the model with a predefined amplitude and time course and was modulated, in the subsequent simulations, as to represent a progressive weakening of the muscle or, at the contrary, an exceeding activation. Ground reaction forces and moments obtained from gait analysis were applied to the foot. To analyze the effects of changing the RF force, we first computed the moments produced by its complementary synergistic muscles: the VA at the knee, and the Iliopsoas at the hip. This was done by running the simulation in an inverse dynamics mode, with normal kinematics and RF force as input. Since the RF force produced a moment at hip and knee, the joint moments computed by the model represented the additional contribution required to other muscles than RF. These moments were named the "residual" moments. When the model was used in a forward dynamics mode, the residual moments plus the moments produced by the RF force reproduced exactly the previous kinematics.

To simulate the effect of different RF forces we just changed the RF force in a predetermined time window and we obtained different kinematics. This was the objective of our study.

In this approach, the complementary muscles were supposed to be unaffected by changes of the RF force, so that this condition was considered: "without compensation". We also considered the possibility of compensation by other muscles to the change of RF force. Of course, it is difficult to imagine what strategy the motor control system can adopt. For the sake of simplicity, we assumed that compensation consisted in reproducing a percentage of the moment originally produced by the RF. For example, 0% compensation means that kinematics will be determined by the residual moments only; 50% compensation means that the residual moments will be incremented by 50% of the moments originally produced by RF; 100% compensation means that the moments originally produced by RF are entirely added to the residual moments so that the resulting kinematics are exactly the original ones.

This criterion of simulating compensations to the altered RF force was applied in two different conditions: hip joint constrained by the predefined, normal kinematics (this represents a full compensation at the hip joint) and both hip and knee free to move under the effect of the forces and moments applied. compensation at the hip joint) and both hip and knee free to move under the effect of the forces and

*Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 4 of 12

According to the purpose of our study, RF force alterations were applied in the time window from 50% to 80% of the gait cycle, which corresponded to the late stance-early swing phase. moments applied. According to the purpose of our study, RF force alterations were applied in the time window from 50% to 80% of the gait cycle, which corresponded to the late stance-early swing phase.

#### (b) Simulations (b) Simulations

At first, we defined a hypothetical RF force for the time window of interest. The force predicted by the Anybody model, based on static optimization, was deemed to be appropriate. The time course was in line with the activation profiles of the RF muscle described in the literature [3,4,15], although the peak amplitude was relatively higher, and in our simulations produced residual moments close to zero or even negative. We then downscaled the RF force as to have a peak of about 300 N corresponding to the peak value reported in Anderson and Pandy 2001 [16]. In this way, the moment produced by RF at the knee was about 40% of the total knee extensor moment (Figure 1), while the residual 60% was produced by the complementary muscles (the VA). At first, we defined a hypothetical RF force for the time window of interest. The force predicted by the Anybody model, based on static optimization, was deemed to be appropriate. The time course was in line with the activation profiles of the RF muscle described in the literature [3,4,15], although the peak amplitude was relatively higher, and in our simulations produced residual moments close to zero or even negative. We then downscaled the RF force as to have a peak of about 300 N corresponding to the peak value reported in Anderson and Pandy 2001 [16]. In this way, the moment produced by RF at the knee was about 40% of the total knee extensor moment (Figure 1), while the residual 60% was produced by the complementary muscles (the VA).

**Figure 1.** Above: different percentages of rectus femoris (RF) force input into the model (the force exerted in the unperturbed stride movement was 100%); middle and low: Knee and hip joint moments produced in unperturbed condition by all muscles (Kmom, Hmom, blue line), by the RF muscle (RFmom), and by other muscles than the RF (residual moments). **Figure 1.** Above: different percentages of rectus femoris (RF) force input into the model (the force exerted in the unperturbed stride movement was 100%); middle and low: Knee and hip joint moments produced in unperturbed condition by all muscles (Kmom, Hmom, blue line), by the RF muscle (RFmom), and by other muscles than the RF (residual moments).

Several simulations were then carried out by imposing the unperturbed kinematics at the hip and letting the knee rotate freely. When the residual knee moment and the RF force (100%) were applied all together, the resulting kinematics were superimposable to the one obtained with imposed kinematics (unperturbed stride movement). Next, in a time window from 50% to 80% of the stride time, we set the RF force to zero in order to simulate a complete elimination of the RF muscle, without Several simulations were then carried out by imposing the unperturbed kinematics at the hip and letting the knee rotate freely. When the residual knee moment and the RF force (100%) were applied all together, the resulting kinematics were superimposable to the one obtained with imposed kinematics (unperturbed stride movement). Next, in a time window from 50% to 80% of the stride time, we set the RF force to zero in order to simulate a complete elimination of the RF muscle, without compensation

from any other knee extensor muscle. Then, we simulated a partial compensation by adding an RF force corresponding to 50%, 75%, 90%, and 95% of the original RF force. To better understand the role of the RF, we also simulated a hyper-activity of this muscle by setting its force to 105%, 110%, 125%, and 150% of the original force. The same sequence of simulations was then repeated by letting both the knee and the hip joints free to rotate. (a) Inverse dynamics: The RF moment at the knee joint obtained from the Anybody computation was in agreement with the EMG activity known from other studies [4,6,8]. The power of RF muscle was also computed and is reported in Figure 2 (red line). This variable synthetically represents the energetic contribution of RF to the knee joint movement. Figure 2 shows that in a normal situation the power is negative

the shank and was involved in active knee extension. The power of VA was computed in two different

*Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 5 of 12

understand the role of the RF, we also simulated a hyper-activity of this muscle by setting its force to 105%, 110%, 125%, and 150% of the original force. The same sequence of simulations was then

repeated by letting both the knee and the hip joints free to rotate.

#### **3. Results** (energy is absorbed) from 40 to 71% of the gait cycle in correspondence with the knee flexion occurring at the late stance-early swing phase. From 75 to 85% of the gait cycle where the extension

(a) Inverse dynamics: of the knee started, the power became positive. This indicates that the RF was transferring energy to

The RF moment at the knee joint obtained from the Anybody computation was in agreement with the EMG activity known from other studies [4,6,8]. The power of RF muscle was also computed and is reported in Figure 2 (red line). This variable synthetically represents the energetic contribution of RF to the knee joint movement. Figure 2 shows that in a normal situation the power is negative (energy is absorbed) from 40 to 71% of the gait cycle in correspondence with the knee flexion occurring at the late stance-early swing phase. From 75 to 85% of the gait cycle where the extension of the knee started, the power became positive. This indicates that the RF was transferring energy to the shank and was involved in active knee extension. The power of VA was computed in two different conditions: when RF was contributing as a part of the quadriceps muscle (green line in Figure 2) and without the contribution of RF (paralyzed muscle, blue line). The sequence of negative and positive power observed at early stance phase (out of interest in this study) was dominated by the VA contribution, so that no difference appears between presence or absence of RF (the green band is completely hidden by the blue band). In the phase of interest (late stance-early swing) the negative power of the VA was very small when RF was present (green line) and became higher, in absolute, than the power required by the RF in normal conditions when RF was paralyzed (blue line). conditions: when RF was contributing as a part of the quadriceps muscle (green line in Figure 2) and without the contribution of RF (paralyzed muscle, blue line). The sequence of negative and positive power observed at early stance phase (out of interest in this study) was dominated by the VA contribution, so that no difference appears between presence or absence of RF (the green band is completely hidden by the blue band). In the phase of interest (late stance-early swing) the negative power of the VA was very small when RF was present (green line) and became higher, in absolute, than the power required by the RF in normal conditions when RF was paralyzed (blue line). The missing negative power of the paralyzed RF was replaced by a greater negative power of the VA. Between 75 and 87% of the gait cycle. Further, there was a positive power at the knee joint produced by the VA, which replaced the activity of the paralyzed RF. This indicates that the positive muscle power of the knee extensors was necessary during this phase. It underlines that the initial extension of the knee in swing phase at 75% of the gait cycle (which corresponds to muscular power generation) was not only a passive motion induced by a double pendulum mechanism but also depended on knee extensor action, mainly the one of RF in normal conditions, or of VA in case of a paralyzed RF.

**Figure 2.** Inverse dynamic simulation. Red line = power of RF with normal activity. Green line = power of all vasti with normal RF activity. Blue line = power of all vasti when RF is paralyzed (note that green line and blue line overlap during the stance phase, since RF power is negligible in this phase, and the vasti do not change their activity. Encircled: the need for positive power at the beginning of knee extension is evidenced. Below: Knee joint angle along the stride cycle, where flexion is positive and extension is negative. **Figure 2.** Inverse dynamic simulation. Red line = power of RF with normal activity. Green line = power of all vasti with normal RF activity. Blue line = power of all vasti when RF is paralyzed (note that green line and blue line overlap during the stance phase, since RF power is negligible in this phase, and the vasti do not change their activity. Encircled: the need for positive power at the beginning of knee extension is evidenced. Below: Knee joint angle along the stride cycle, where flexion is positive and extension is negative.

The missing negative power of the paralyzed RF was replaced by a greater negative power of the VA. Between 75 and 87% of the gait cycle. Further, there was a positive power at the knee joint produced by the VA, which replaced the activity of the paralyzed RF. This indicates that the positive muscle power of the knee extensors was necessary during this phase. It underlines that the initial extension of the knee in swing phase at 75% of the gait cycle (which corresponds to muscular power generation) was not only a passive motion induced by a double pendulum mechanism but also depended on knee extensor action, mainly the one of RF in normal conditions, or of VA in case of a paralyzed RF. (b) Forward modeling We first considered the condition of hip kinematics unchanged (full compensation at the hip by the hip muscles, Figure 3). This of course is a hypothetical situation that only can be reproduced in

simulation, but it helps understanding the separate effect of RF at the knee joint. With the RF in place

*Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 6 of 12

#### (b) Forward modeling (100% RF force and residual knee moments applied), the model performed a normal step as expected.

We first considered the condition of hip kinematics unchanged (full compensation at the hip by the hip muscles, Figure 3). This of course is a hypothetical situation that only can be reproduced in simulation, but it helps understanding the separate effect of RF at the knee joint. With the RF in place (100% RF force and residual knee moments applied), the model performed a normal step as expected. After removal of the RF and without any compensation from other knee extensor muscles (0% RF force and residual moments applied), exaggerated knee flexion and foot rising occurred in swing, which prevented the foot from getting in contact with the ground at the end of the stride cycle (Figure 3, left panel). The same occurred at 50% of RF force. At 75% of RF force, the step was still incomplete. In clinical terms this would correspond to abnormal foot contact occurring with the toes (initial toe contact). At 90% and 95% of the RF force, the foot recovered a position, at the end of the stride cycle, that was compatible with foot plant initial contact and heel initial contact, respectively. After removal of the RF and without any compensation from other knee extensor muscles (0% RF force and residual moments applied), exaggerated knee flexion and foot rising occurred in swing, which prevented the foot from getting in contact with the ground at the end of the stride cycle (Figure 3, left panel). The same occurred at 50% of RF force. At 75% of RF force, the step was still incomplete. In clinical terms this would correspond to abnormal foot contact occurring with the toes (initial toe contact). At 90% and 95% of the RF force, the foot recovered a position, at the end of the stride cycle, that was compatible with foot plant initial contact and heel initial contact, respectively. When the RF was hyperactive, a reduced knee flexion was observed (Figure 3, right panel) that in some cases prevented the foot clearance (the forefoot moved below the ground level). As shown in Figure 3, knee extension was faster than at basic condition, so that at the time corresponding to stride duration the knee was hyperextended and the forward displacement of the foot was exaggeratedly enhanced.

**Figure 3.** Left panel: kinematics resulting from decreasing RF force; right panel: kinematics resulting from increasing RF force. The model is represented in the position achieved at the end of stride time (100% stride time). The red and blue curves represent the trajectory of a point located on the heel of, respectively, the right side (unaffected, blue line) and the left side (our test limb, red line). The hip joint was fully compensated (hip kinematics imposed). **Figure 3.** Left panel: kinematics resulting from decreasing RF force; right panel: kinematics resulting from increasing RF force. The model is represented in the position achieved at the end of stride time (100% stride time). The red and blue curves represent the trajectory of a point located on the heel of, respectively, the right side (unaffected, blue line) and the left side (our test limb, red line). The hip joint was fully compensated (hip kinematics imposed).

When the RF was hyperactive, a reduced knee flexion was observed (Figure 3, right panel) that in some cases prevented the foot clearance (the forefoot moved below the ground level). As shown in Figure 3, knee extension was faster than at basic condition, so that at the time corresponding to stride duration the knee was hyperextended and the forward displacement of the foot was exaggeratedly enhanced.

Figure 4 depicts the knee flexion angles for all the simulated percentages of RF force compensation (50%, 75%, 90%, 95%, 100%) and RF hyperactivation (105%, 110% 125%, 150%). It clearly reveals that knee flexion increases as far as the RF force decreases, and, at the contrary, knee flexion decreases while the RF force increases. *Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 7 of 12 clearly reveals that knee flexion increases as far as the RF force decreases, and, at the contrary, knee flexion decreases while the RF force increases.

**Figure 4.** Changes of the knee joint angle with increasing and decreasing RF force. The hip joint was fully compensated (hip kinematics imposed). **Figure 4.** Changes of the knee joint angle with increasing and decreasing RF force. The hip joint was fully compensated (hip kinematics imposed).

Similar effects of RF removal and hyper-activation were obtained when the hip joint as well as the knee joint was left free to rotate under the effect of the RF force and residual hip moment (hip kinematics were not imposed). The entire lower limb kinematics dramatically changed for large reductions in RF force (condition 50% and 75% in Figure 5). Similar effects of RF removal and hyper-activation were obtained when the hip joint as well as the knee joint was left free to rotate under the effect of the RF force and residual hip moment (hip kinematics were not imposed). The entire lower limb kinematics dramatically changed for large reductions in RF force (condition 50% and 75% in Figure 5).

Less apparent deviations were observed for RF force enhancement. For example, the increase in stride length seemed not too big at 150% RF (Figure 5, right panel). However, it has to be noted that the foot largely moves below ground level, and this is related to reduced movement at both the hip and knee. In clinical terms, this would correspond to great difficulty to obtain foot clearance. Figure 6 shows the changes in hip and knee joint angles. Decreasing the RF force produced an increase in both hip and knee flexion; the opposite occurred with increasing RF force. Compared to Figure 4, the reduction in knee flexion was more pronounced in this case in which the RF also had an effect at the hip.

**Figure 4.** Changes of the knee joint angle with increasing and decreasing RF force. The hip joint was

Similar effects of RF removal and hyper-activation were obtained when the hip joint as well as

*Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 7 of 12

clearly reveals that knee flexion increases as far as the RF force decreases, and, at the contrary, knee

flexion decreases while the RF force increases.

fully compensated (hip kinematics imposed).

reductions in RF force (condition 50% and 75% in Figure 5).

**Figure 5.** Left: kinematics resulting from decreasing RF force; right panel: kinematics resulting from increasing RF force. The model is represented in the position achieved at the end of stride time (100% stride time). The red and blue curves represent the trajectory of a point located on the heel of, respectively, the right side (unaffected, blue line) and the left side (our test limb, red line). The hip joint flexion was not imposed (no compensation for RF force changes). and knee. In clinical terms, this would correspond to great difficulty to obtain foot clearance. Figure 6 shows the changes in hip and knee joint angles. Decreasing the RF force produced an increase in both hip and knee flexion; the opposite occurred with increasing RF force. Compared to Figure 4, the reduction in knee flexion was more pronounced in this case in which the RF also had an effect at the hip.

**Figure 6.** Changes of the hip and knee joint angles with increasing and decreasing RF force. The kinematics of the hip and knee were not imposed and resulted from the application of the RF force plus the residual moments (see text). **Figure 6.** Changes of the hip and knee joint angles with increasing and decreasing RF force. The kinematics of the hip and knee were not imposed and resulted from the application of the RF force plus the residual moments (see text).

simulations at the time where the foot was in a proper position to contact the ground. Figure 7 shows

the respective foot positions, and Figure 8 reports the step parameters.

Figures 3 and 5 show that when the RF force is decreased, more time would be required for the

Figures 3 and 5 show that when the RF force is decreased, more time would be required for the foot to reach the ground. Therefore, we allowed the simulation to continue beyond the cycle time without any internal joint moment in order to get the position of the foot (angle and step length) at the foot–floor contact time. We also computed the same parameters in the "increased RF force" simulations at the time where the foot was in a proper position to contact the ground. Figure 7 shows the respective foot positions, and Figure 8 reports the step parameters. *Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 9 of 12 *Appl. Sci.* **2020**, *10*, x FOR PEER REVIEW 9 of 12

**Figure 7.** Foot position achieved at foot-floor contact when RF force was changed (the numbers show the percentage of RF force). The two columns on the left refer to a condition in which only the knee joint was controlled by the RF force, while the hip joint kinematics was imposed (this corresponded to full compensation by hip muscles); the two columns on the right refer to a condition in which both hip and knee were dependent on RF force. **Figure 7.** Foot position achieved at foot-floor contact when RF force was changed (the numbers show the percentage of RF force). The two columns on the left refer to a condition in which only the knee joint was controlled by the RF force, while the hip joint kinematics was imposed (this corresponded to full compensation by hip muscles); the two columns on the right refer to a condition in which both hip and knee were dependent on RF force. **Figure 7.** Foot position achieved at foot-floor contact when RF force was changed (the numbers show the percentage of RF force). The two columns on the left refer to a condition in which only the knee joint was controlled by the RF force, while the hip joint kinematics was imposed (this corresponded to full compensation by hip muscles); the two columns on the right refer to a condition in which both hip and knee were dependent on RF force.

**Figure 8.** Changes of step parameters when the RF force was changed (the percentage of RF force is reported on the abscissa). Left panel refers to conditions in which the hip kinematics was imposed; right panel: both the hip and knee kinematics were determined by RF force. From top to bottom: time at initial contact; step length; height of the heel at time 100% of stride cycle; angle of the foot plant in **Figure 8.** Changes of step parameters when the RF force was changed (the percentage of RF force is reported on the abscissa). Left panel refers to conditions in which the hip kinematics was imposed; right panel: both the hip and knee kinematics were determined by RF force. From top to bottom: time at initial contact; step length; height of the heel at time 100% of stride cycle; angle of the foot plant in relation to the horizontal line (negative angles mean toes down). **Figure 8.** Changes of step parameters when the RF force was changed (the percentage of RF force is reported on the abscissa). Left panel refers to conditions in which the hip kinematics was imposed; right panel: both the hip and knee kinematics were determined by RF force. From top to bottom: time at initial contact; step length; height of the heel at time 100% of stride cycle; angle of the foot plant in relation to the horizontal line (negative angles mean toes down).

force above 120%, when both hip and knee were free to move, because they stayed almost fixed, and

relation to the horizontal line (negative angles mean toes down).

the foot soon penetrated the ground.

the foot soon penetrated the ground.

Note that these figures do not report the results of when the RF force reduced to 50%, since the foot did not reach the ground. As well, it was not possible to calculate the step parameters for the RF force above 120%, when both hip and knee were free to move, because they stayed almost fixed, and the foot soon penetrated the ground.

#### **4. Discussion**

The role of the RF as a bi-articular muscle during the gait cycle even in normal human walking is difficult to understand. The effect of this muscle on the mechanics of swing is even more difficult to estimate as the segments of the leg are accelerated in pre- and early swing while muscle activity during the later swing phase is low. The anatomical situation suggesting that the RF acts as a hip flexor and a knee extensor at the same time, and also the possibility that other muscles can compensate for its weakness [17] makes it difficult to make conclusions about the effects of this muscle on each joint. Several studies based on EMG [6–8,11] and musculoskeletal modelling [1,2,18,19] have revealed that the RF is important during the transition from stance phase to swing phase to control knee flexion. Its detailed function during this period and its consequence on gait kinematics especially during swing phase were the main focus of this study.

At the knee level, the RF indeed controlled knee flexion from 40 to 75% of the gait cycle. This was achieved, as shown in our inverse dynamics analysis, by absorbing mechanical power and breaking the knee flexion movement. Without knee extensors, the knee would flex more extensively. After the peak of knee flexion, however, RF started to produce a positive power, which resulted in knee extension from 75 to 85% of the gait cycle. After this push, no more muscle power was necessary to further extend the knee because the acceleration of the leg followed a double pendulum movement controlled by negative hamstring power. The importance of the initial push for knee extension during swing was described for stiff knee gait but without considering the effect on foot position at terminal swing [10].

The forward simulation indeed showed that the RF accelerated the shank segment into an extension movement of the knee. Thus, the RF controls step length and is one factor that determines an adequate foot prepositioning. This has clinical significance, as with reduced knee extensor force (RF removal or weakening), the knee remains flexed at late swing, which results in a short step and initial toe contact. The simultaneous effect of the RF on the hip, however, had little influence when the activity was reduced. In contrast, hyperactivity led to more knee stiffness if the RF had effects on both joint levels, leading to a foot clearance problem. These forward modelling results are of clinical interest, even though our model, of course, can only simulate the mechanical aspects and not the control strategies that a patient would adopt to compensate the RF alterations. We have shown that stiffness at the knee during the swing phase, at least in a healthy person, depends on the increased RF action at both joint levels (although more strongly at the knee). As RF activity during gait has been previously described [3–8], we limited our study to its mechanical effect on swing phase kinematics in principle, which required gait data from one normal person only as an example. The amount of the values may change with changing segment sizes and proportions. However, the anatomy is similar and the gait pattern changes little after the age of four [12]. The effect of bi-articular muscles during gait depends on the anatomy but also on the interplay with other muscles crossing the joints and on physical parameters such as acceleration and inertia. In a clinical situation with pathological gait, the situation is even more complicated. We demonstrated the effects of lower and higher than normal RF activities on gait mechanics, especially with reference to the foot position. RF accelerates the leg segments in late stance-early swing and hence it can affect the whole subsequent movements: an adequate knee extension at late swing and preparation of the leg for the initial heel contact. As shown in our simulations, RF weakness may result in toe walking. This effect may be considered especially for patients who already have problems with knee extension during swing phase and exhibit initial toe contact. Hyperactivity, in contrast, can hinder foot clearance. Both situations occur in clinical practice. The new knowledge may help to identify RF as a possible factor. More sophisticated diagnostic tools

are required to facilitate the process of treatment decision-making. Incorporation of musculoskeletal modelling into clinical work may be one promising solution.

#### **5. Conclusions**

In conclusion, our study confirms the important role of RF in controlling the swing phase of gait. Its activity must be carefully controlled during late stance-early swing to prevent abnormal knee flexion and improper foot positioning at the end of the swing phase. Reducing the RF contribution may produce an increased knee flexion and an inappropriate foot position at foot–ground contact; increased RF activity reduces knee flexion and hinders foot clearance. Like with every forward modelling study, the results may be less apparent in a clinical situation as motor control can produce a compensatory muscle activation pattern. Nevertheless, our results on the mechanical effects of the enhanced or reduced RF strength showed detrimental effects on knee joint excursion, foot placement, and spatio-temporal parameters.

**Author Contributions:** All authors contributed to this paper: C.A.F. performed the forward modelling computations and analysis, C.W. the inverse dynamics part on the gait data, R.B. was responsible for the research question and supervision of the manuscript preparation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Development of a Musculoskeletal Model of Hyolaryngeal Elements for Understanding Pharyngeal Swallowing Mechanics**

**Takuya Hashimoto 1,\* , Mariko Urabe <sup>1</sup> , Foo Chee-Sheng <sup>1</sup> , Atsuko Murakoshi <sup>2</sup> , Takahiro Kikuchi <sup>3</sup> , Yukihiro Michiwaki <sup>3</sup> and Takuji Koike <sup>2</sup>**


Received: 14 August 2020; Accepted: 7 September 2020; Published: 9 September 2020

**Featured Application: The potential application of this study is to provide a key tool to help researchers and clinicians understand how neuro, muscular, and skeletal systems are involved in swallowing. It would contribute to the clinical decision-making process in treating dysphagia associated with neuromuscular disease.**

**Abstract:** A detailed understanding of muscle activity in human swallowing would provide insights into the complex neuromuscular coordination underlying swallowing. The purpose of this study was to introduce musculoskeletal analysis to investigate muscle activities involved in swallowing as there are limitations on studying comprehensive muscle activation patterns by conventional methods such as electromyography (EMG) measurement. A musculoskeletal model of swallowing was newly developed based on the skeletal model made from CT data of a healthy volunteer. Individual muscle forces were predicted in pharyngeal swallowing by inverse dynamics' computations with static optimization, in which the typical trajectories of the hyoid bone and thyroid cartilage analyzed from videofluoroscopic (VF) data of the volunteer were used. The results identified the contribution of individual muscles in pharyngeal swallowing in relation to the movements of the hyoid bone and thyroid cartilage. The predicted sequence of muscle activity showed a qualitative agreement with salient features in previous studies with fine wire EMG measurements. This method, if validated further by imaging and EMG studies, enables studying a broader range of neuromuscular coordination in swallowing. The proposed method offers an avenue to understanding the physiological mechanisms of swallowing and could become useful to evaluate rehabilitation effects on dysphagia.

**Keywords:** deglutition; musculoskeletal model; muscle activity; simulation; swallowing biomechanics

#### **1. Introduction**

The increase in the number of aged in the Japanese population is faster than any other country and the elderly (aged 65 or over) comprised 27.3% in 2016 [1]. One critical issue accompanying aging is impairment of physical and cognitive abilities, with the decline in the ability to swallow particularly causing serious problems including aspiration pneumonia. Therefore, dysphagia difficulty with swallowing is a serious health concern in an aging society, and pneumonia has been the third most common cause of death in Japan from 2011 [2], with most of the pneumonia in the elderly ascribed to

aspiration pneumonia caused by dysphagia [3,4]. Disability in swallowing also increases the risks of suffocation, undernutrition, and dehydration, which have the potential to decrease the quality of life. Therefore, prevention and treatment of dysphagia are a social issue in an aging society with extended healthy life spans and social welfare concerns.

Swallowing, or deglutition, is a coordinated and smoothly functioning process to propel food and drink from the mouth to the esophagus. Since it involves complex neuromuscular coordination [5], it is difficult to understand the physiological mechanisms. Videofluoroscopic (VF) examination [6,7] is commonly used in research and clinical settings as the preferred standard method for assessing swallowing function because of the comprehensive information it provides. This method allows anatomic and kinematic analysis such as the movement of food bolus and swallowing-related organs. Videoendoscopic (VE) examination [8,9] is also widely used as a reliable and portable technique to observe the pharyngolaryngeal complex directly. However, the information on muscle activation patterns during swallowing cannot be obtained by these methods, though swallowing is controlled by the coordination of many muscles [5]. Surface electromyograms (sEMG) [10,11] are the best way to estimate muscle activation patterns non-invasively so far, but the simultaneous measurement of swallowing-related muscles is limited. The main reason for this is that most swallowing-related muscles are located deep below the surface and overlap as if forming a layered structure. As a result, sEMG cannot avoid cross-talk problems in the simultaneous measurement of many muscles in swallowing. Therefore, conventional methods cannot provide information to quantify the sequential pattern of the muscles in swallowing.

It is impossible to measure internal forces, such as joint contact forces and muscle strength, non-invasively with conventional techniques. Hence, biomechanical analysis using musculoskeletal models [12–15] has been considered as an alternative technique to estimate such internal forces. For example, musculoskeletal models for the lower extremities were used to analyze complicated inter-muscle coordination in dynamic movement, such as in walking [16,17], running [18], and jumping [19]. Upper extremity models [20,21] were also developed for the biomechanical analysis of dynamic motion of the upper limbs, such as in throwing [22], propelling of wheelchair [23], and so forth. Presently, large-scale musculoskeletal analysis is readily conducted, as software for developing a complicated musculoskeletal model has become available commercially [12,13] or as open-source software [14]. It may be expected that muscle activity during swallowing can be estimated in a similar manner as dynamic body movements even though most of the muscle activity involved in swallowing cannot be determined in vivo. However, there are few studies reporting musculoskeletal models of swallowing [24,25]. The authors have reported a musculoskeletal model of swallowing by analyzing normal swallowing motion considering muscle activity [26,27], but an analysis of details of the muscle activity in swallowing is still difficult and challenging.

Previous studies on computational biomechanical analysis of swallowing mainly focused on the interaction between the bolus and organs considering fluid dynamics [28–30]. In these studies, the deformation of organs was preliminarily provided by measured results and the activity of the muscles was not considered. However, estimating muscle activity would contribute to an understanding of the physiological aspects of swallowing as muscle expansion and contraction affect the deformation of organs.

This paper presents the design of a musculoskeletal model of swallowing and an investigation of muscle activity during pharyngeal swallowing. The skeletal model considers the hyoid bone and thyroid cartilage as a central element in the study of the swallowing function. The configuration of the swallowing-related muscles is defined based on anatomical knowledge, and the muscles are represented as a Hill-type model. The trajectories of the hyoid bone and thyroid cartilage using videofluoroscopic data of a healthy subject were analyzed for the dynamics' computations. The forces of individual muscles were estimated through inverse dynamics' computations with static optimization to study the activation pattern of the muscles and their contributions in the swallowing. The musculoskeletal

analysis is one method that can be used in the study of details of neuromuscular coordination of swallowing in future work.

#### **2. Materials and Methods**

#### *2.1. Motion Analysis of the Hyoid Bone and Thyroid Cartilage during Swallowing*

The normal swallowing process is generally classified into three phases according to the location of bolus: an oral phase, a pharyngeal phase, and an esophageal phase. The pharyngeal phase, where the bolus is transported from the pharynx into the esophagus, is especially important. Here, the larynx and pharynx move in a superior direction and the tongue base moves in a postero-superior direction in concert with movements of the hyoid bone and thyroid cartilage. As a result, the epiglottis folds down to prevent the bolus from entering the trachea leading it to the esophagus [31]. If this process is not completed fully, food and drink will enter the trachea and cause aspiration. Analyzing the movements of the hyoid bone and thyroid cartilage during the pharyngeal phase is important in the evaluation of swallowing because the movement of the larynx is affected by muscles attached to the hyoid bone and thyroid cartilage.

#### 2.1.1. Method of the Analysis

To analyze two-dimensional movements of the hyoid bone and thyroid cartilage, we used an X-ray video (DR-2000F, Hitachi Medical Co., Tokyo, Japan) obtained from a videofluoroscopic (VF) examination, which recorded normal swallowing with a healthy, 25-year-old male who was not diagnosed with any other eating disorder or neurological diseases. Here, the volunteer swallowed 5 mL of water, including the contrast of 2.21 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>2/s viscosity, once. The volunteer gave his consent before the VF examination. The experimental procedure was approved by the ethics committee of Japanese Red Cross Musashino Hospital (No. 705). In the analysis, as shown in Figure 1, we manually extracted the region of the hyoid bone and the vocal ligament as the part of thyroid cartilage, frame by frame. The displacement and posture of each skeletal element were measured in relation to the initial position and orientation. Then, the measured displacements and postures were transformed to a global coordinate system with the center at the left–right occipital condyles. Each skeletal element has a local coordinate system at its center of gravity, as shown in Figure 1, and the volume of the two elements was estimated from skeletal data of a CT scan (Revolution CT, GE Healthcare Japan Co., Tokyo, Japan) of the same male. In the study here, the cricoid cartilage was included as part of thyroid cartilage. 2.21 × 10ି m<sup>ଶ</sup>⁄s

**Figure 1.** Trajectory analysis of the hyoid bone and thyroid cartilage. (**a**) Analyzed videofluoroscopic (VF) image, (**b**) global and local coordinate systems.

#### 2.1.2. Results of the Analysis

<sup>ଵ</sup> <sup>ଶ</sup>

s

<sup>ୱ</sup> ୣ ୱ s <sup>ଵ</sup> s Figure 2 shows the trajectories of the hyoid bone and thyroid cartilage in the sagittal plane, smoothed by a simple moving average. Figure 3 shows the displacements of the hyoid bone and thyroid cartilage in the *X*- and *Y*-axis directions, where *P*<sup>s</sup> and *P*<sup>e</sup> represent the starting and ending points, respectively. These figures show that the movements of the two elements can be divided into three phases. In the first phase, between *P*<sup>s</sup> (0.90 s) and *P*<sup>1</sup> (1.36 s), the later oral phase, the hyoid bone

moves upward and slightly backward while thyroid cartilage moves upward and slightly forward. Here, the bolus is held in the oral cavity and its head reaches the soft palate. In the second phase, between *P*<sup>1</sup> and *P*<sup>2</sup> (1.73 s), the earlier pharyngeal phase, the two move together in an antero-superior direction and the soft palate is elevated. The bolus enters the oropharynx and the closure of the larynx by the epiglottis occurs. Then, the bolus moves to the hypopharynx and the esophagus. In the last phase, between *P*<sup>2</sup> and *P*<sup>e</sup> (2.23 s), the later pharyngeal phase, they move in a postero-inferior direction to return to the resting positions, slightly different from the initial positions. The movement of the hyoid bone exhibits a typical triangular trajectory, as described in the literature [32]. <sup>ଶ</sup> ୣ s <sup>ଶ</sup> ୣ s

**Figure 2.** Trajectories of the hyoid bone and thyroid cartilage. (**a**) Hyoid bone, (**b**) thyroid cartilage.

**Figure 3.** Displacements of the hyoid bone and thyroid cartilage. (**a**) Displacement in *X*-axis, (**b**) displacement in *Y*-axis, (**c**) rotation angle around *Z*-axis.

#### *2.2. Development of the Musculoskeletal Model of Swallowing*

#### 2.2.1. Modeling of the Swallowing-Related Muscles

To analyze the muscle activity during swallowing, we developed a musculoskeletal model where a muscle is represented as an active constrictive wire actuator, by reference to the previous study [15]. In this model, the following 12 kinds of muscles were included as muscles acting on the hyoid bone and thyroid cartilage, as shown in Figures 4 and 5: the geniohyoid muscle, the stylohyoid muscle, the digastric muscle, the mylohyoid muscle, the sternohyoid muscle, the omohyoid muscle, the middle constrictor muscle, the thyrohyoid muscle, the sternothyroid muscle, the inferior constrictor muscle, the cricopharyngeus muscle, and the stylopharyngeus muscle. The muscles on the left and right sides are distributed symmetrically about the sagittal plane for simplicity. The origin and insertion of the muscles are shown in Table 1, defined based on anatomical description [33] and CT data in addition to BodyParts3D [34,35]. Broad muscles such as the mylohyoid and pharyngeal constrictor muscles are represented as multiple wires. Muscles with two muscle bellies such as the digastric and omohyoid muscles are modeled as two separate muscles with an intermediate tendon. The modeling of the musculoskeletal model was done using 3D modeling software, 3ds Max®. Further detailed description of the modeling of the swallowing-related muscles is as documented in previous studies [26,27].

**Figure 4.** Swallowing-related muscles with wire models. Left image is an anatomical image of the muscle, and the wire model is shown on the right. (**a**) Geniohyoid (GH), (**b**) stylohyoid (SH) and digastric (AD, PD), (**c**) mylohyoid (AMH, PMH), (**d**) sternohyoid (SN) and omohyoid (OH), (**e**) middle pharyngeal constrictor (PG1-PG4, PL1-PL4), (**f**) thyrohyoid (TH) and sternothyroid (ST), (**g**) inferior pharyngeal constrictor (PT1-PT4) and cricopharyngeus (CP1, CP2), (**h**) stylopharyngeus (SP).

**Figure 5.** Overview of the musculoskeletal model of swallowing. The left image is a frontal view and the right image is a lateral view. The black lines represent muscle models and the blue lines represent membrane springs.

#### 2.2.2. Modeling of Damping and Stiffness of Biological Tissue

× ି × ି To incorporate the viscous-elastic properties of biological tissue, such as skin, fat, muscle, membranes, and others, imaginary linear viscous and elastic components were added to constrain the movements of the hyoid bone and thyroid cartilage in the *X*- and *Y*-axis directions. However, it was quite difficult to represent the net effect accurately because the measurement of the complex mechanical properties of human soft tissue is not straightforward. Therefore, we assumed that the elastic constant of soft tissue for the movements in the *X*- and *Y*-axis directions was *k* = 25.0 N/m based on the previous study [36]. Then, the viscosity constant *c* was calculated to satisfy overdamping condition (*c* > √ *mk*) due to high damping effect of the soft tissue. Here, *m* was the mass of the hyoid bone or thyroid cartilage. As a result, we assumed *c* = 2.0 Ns/m.

Further, the thyrohyoid membrane between the hyoid bone and thyroid cartilage, as well as the membrane between the cricoid cartilage and trachea, were represented as multiple linear elastic components. The elastic constant of each membrane component was experimentally defined as *k<sup>m</sup>* = 20.0 N/m through an iterative process, so as to converge quadratic optimization program described below.


**Table 1.** List of swallowing-related muscles and their physiological cross-sectional areas (PCSAs).

PCSA: Physiological cross-sectional area (\*) represents the estimated value, which is not validated.

#### 2.2.3. Modeling of the Mechanical Properties of Muscles

In swallowing, the hyoid bone and thyroid cartilage are mainly moved by a tug-of-war among the muscles, and the passive force of a muscle generated when it is stretched in addition to the active contractile force must be considered. The Hill-type muscle model [37], which is commonly used in musculoskeletal studies, was used as a simple phenomenological computational model to represent the mechanical properties of the muscles in this study.

The Hill-type muscle model consists of a contractile element (CE), a parallel spring element (PE), and a serial spring element (SE). The capacity to generate the force of a muscle is defined by the force-length and force-velocity properties of the contractile element and the nonlinear spring properties of the passive elements. The total muscle force is the sum of the active (*Fce*) and passive forces (*Fpe*) as below:

$$F\_{mus} = F\_{c\varepsilon} + F\_{p\varepsilon}.\tag{1}$$

As reported elsewhere [38,39], the active force generated by the contractile elements can be described by the following Equations (2) and (3):

$$F\_{\rm ce}(a, l\_{\rm ce}, \dot{l}\_{\rm ce}) = a f\_{\rm lce}(l\_{\rm ce}) f\_{\rm lce}(\dot{l}\_{\rm ce}) \cos(a) f\_{\rm max} \tag{2}$$

$$f\_{\text{max}} = \sigma\_{\text{max}} \text{PCSA} \tag{3}$$

where *a* is the level of activation of a muscle ranging from 0 and 1, *flce*(*lce*) represents the relationship between the normalized muscle force and its length, *<sup>f</sup>vce* . *lce* represents the relationship between the normalized muscle force and its contractile velocity, α is the pennation angle of the muscle, σ*max* is the maximum isometric muscle stress, and *PCSA* is the physiological cross-sectional area of the muscle. The force–length relationship, *flce*(*lce*), shows a Gaussian distribution around the optimal muscle fiber length *lce*<sup>0</sup> as below:

$$f\_{lcc}(l\_{cc}) = \exp\left[-\left(\frac{l\_{cc} - l\_{cc0}}{l\_{csh}}\right)^2\right].\tag{4}$$

Here, *lce* and *lce*<sup>0</sup> are the current and optimal lengths of a contractile element and *lcesh* is a parameter defining the shape of the Gaussian force–length curve. Note that the optimal length *lce*<sup>0</sup> was assumed as the muscle length in the resting state in this study. The force–velocity relationship, *<sup>f</sup>vce* . *lce* , is formulated as follows: .

$$f\_{\rm cr}(\dot{l}\_{\rm cr}) = \begin{cases} 0 & \left(\dot{l}\_{\rm cr} \le -v\_{\rm max}\right) \\ \frac{V\_{sh}\Big{(v\_{\rm max}(a,l\_{\rm cr}) + \dot{l}\_{\rm cr})}}{V\_{sh}v\_{\rm max}(a,l\_{\rm cr}) - \dot{l}\_{\rm cr}} & \left(-v\_{\rm max} \le \dot{l}\_{\rm cr} < 0\right) \\ \frac{V\_{sh}V\_{sh}v\_{\rm max}(a,l\_{\rm cr}) + V\_{ml}l\_{\rm cr}}{V\_{sh}V\_{sh}v\_{\rm max}(a,l\_{\rm cr}) + \dot{l}\_{\rm cr}} & \left(\dot{l}\_{\rm cr} \le 0\right) \end{cases} \tag{5}$$

$$v\_{\rm max}(a,l\_{\rm cr}) = V\_{\rm vm}\{1 - V\_{\rm cr}[1 - af(l\_{\rm cr})]\} \tag{6}$$

where *vmax* is the maximum contractile velocity of the muscle without load, and the muscle contractile velocity slows as the load increases; *Vsh* and *Vshl* are parameters defining the concavity of the Hill curve during shortening and lengthening, respectively; *Vml* and *Vvm* are the maximum velocities during concentric contraction and during maximum isometric contraction, respectively; and *Ver* is the effect of the activation on the maximum velocity.

The passive force–length relationship of a muscle is represented by an exponential function as below:

$$F\_{\rm pe}(l\_{\rm mus}) = \begin{cases} 0 & (l\_{\rm mus} \le l\_{\rm mus0})\\ \frac{F\_{\rm max}}{\varepsilon^{\rm PE\_{\rm shl}}} \left( e^{\frac{l\_{\rm mus0} - l\_{\rm rms0}}{l\_{\rm rms0}} \cdot \frac{PE\_{\rm shl}}{PE\_{\rm ram}}} - 1 \right) & (l\_{\rm mus} > l\_{\rm mus0}) \end{cases} \tag{7}$$

with *PEsh* and *PExm* as the shape and range parameters of the passive force and *lmus* and *lmus*<sup>0</sup> as the current and optimal muscle lengths. The parameter values of the muscles in this study were defined by referring to previous work [38,39], as shown in Table 2.

#### 2.2.4. Physiological Cross-Sectional Area (PCSA)

The physiological cross-sectional area (PCSA) is used to estimate the maximum muscle force because the maximum muscle force is related to PCSA. Table 1 shows the PCSA of swallowing-related muscles included in the musculoskeletal model, which were defined based on the results of physiological studies [40–42]. However, the PCSA of pharyngeal constrictor and cricopharyngeus muscles (PG, PL, PT, and CP) were estimated based on BodyParts3D [34,35] in addition to anatomical images [33] as these values were not investigated in the previous work. Here, we calculated the PCSAs of pharyngeal constrictor muscles based on their 3D polygon data using 3D modeling software, 3ds Max®.


**Table 2.** Parameters of the muscle model.

#### 2.2.5. Muscle Force Estimation Using an Optimization Program

#### (1) Equation of Motion

Using Lagrangian mechanics and setting the generalized coordinates as *q*, the dynamic equations (such as Equation (8)) of the musculoskeletal model in the following form were obtained:

$$\mathbf{M}\ddot{\boldsymbol{\eta}} + \mathbf{C}\dot{\boldsymbol{\eta}} + \mathbf{K}\boldsymbol{\eta} + \mathbf{F}\_{\text{pu}}^{\text{m}} + \mathbf{g}(\boldsymbol{\eta}) = \mathbf{Q} \tag{8}$$

where *M* is the inertia matrix with a mass *M* and a moment of inertia *J* and *C*and *K*are damping and stiffness matrices, respectively. The *F m pa* is the matrix representing the passive force of the muscles. Thematrix *g* represents the effect of gravity. The *q* = h *x h* , *y h* , *z h* , θ *h* , *x th* , *y th* , *z th* , θ *th*i*<sup>T</sup>* stands for positions and postures of the hyoid bone and thyroid cartilage. The *Q* = [*F h x* , *F h y* , *F h z* , *T h* , *F th x* , *F th y* , *F th z* , *T th*i*<sup>T</sup>* is the generalized force, which includes forces and torques acting on the hyoid bone and thyroid cartilage. The *h* and *th* indices stand for parameters related to the hyoid bone and thyroid cartilage, respectively. The mass of the hyoid bone and thyroid cartilage, *M<sup>h</sup>* and *Mth*, were estimated to be 5.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> kg and 1.5 <sup>×</sup> <sup>10</sup>−<sup>2</sup> kg, respectively, by calculating the volumes from the skeleton data. Further, the moment of inertia about the center of gravity of the hyoid bone and thyroid cartilage, *J h* and *J th*, were estimated using the function of 3D CAD software, SoldWorks®, to be 4.12 <sup>×</sup> <sup>10</sup>−<sup>7</sup> kg m<sup>2</sup> and 2.92 <sup>×</sup> <sup>10</sup>−<sup>6</sup> kg m<sup>2</sup> , respectively.

Using the dynamic Equation (8), inverse dynamics analysis was performed to compute the generalized force *Q* incorporating the forces and torque that realize the motion of the hyoid bone and thyroid cartilage, as shown in Figure 2. Here, the differentiation of *q* with time was computed by a finite-difference method. In addition, the effect of gravity was eliminated (*g* = **0**) for simplicity.

#### (2) Muscle Force Estimation

The developed musculoskeletal model has a redundant number of muscles compared with the degrees of freedom of the hyoid bone and thyroid cartilage. Due to this redundancy, estimating the muscle forces from the obtained generalized forces was an ill-posed problem, and to solve it, Crowninshield's cost function *O*(*f*) [16] was used. The function was formed as a quadratic program as below:

$$O(f) = \sum\_{i=1}^{N} \left(\frac{f\_i}{PCSA\_i}\right)^2 \to \min \tag{9}$$

$$\text{is subject to: } \mathcal{F}\_{\mathbf{x}} = \sum\_{i=1}^{N} f\_{\mathbf{i}\mathbf{x}} \; \mathcal{F}\_{\mathbf{y}} = \sum\_{i=1}^{N} f\_{\mathbf{i}\mathbf{y}} \; \mathcal{F}\_{\theta} = \sum\_{i=1}^{N} f\_{\mathbf{i}} r\_{i} \tag{10}$$

$$f\_i^{\min} \le f\_i \le f\_i^{\max} \quad \text{ ( $i = 1, \dots, N$ )}\tag{11}$$

where *PCSA<sup>i</sup>* is the physiological cross-sectional area (PCSA) of the muscle *i*, *f<sup>i</sup>* is the muscle force, *ri* is the moment arm of muscle *i* at the hyoid bone or thyroid cartilage, *f max i* is the maximum force of the muscle *i*, and *f min i* is the minimum. A combination of possible muscle forces that minimize the total muscle force was obtained by optimizing the cost function (9). Function (10) is the equality constraint and possible muscle forces were required to satisfy the forces and the torque computed by inverse dynamics. Possible muscle forces were also constrained to be equal to or above zero by the inequality constraint (11) because muscles generate only contractile force. Here, *f max <sup>i</sup>* was calculated for each muscle by Equation (3) and *f min <sup>i</sup>* was zero (*f min i* = 0). The quadratic program was solved using MATLAB® function.

#### **3. Results**

#### *3.1. Results of Muscle Force Estimation*

Figure 6 shows the muscle forces (%*Fmax*) estimated by the inverse dynamics computation of swallowing and the static optimization, which were normalized by the maximum forces as below:

$$\%F\_{\text{max}} = \frac{f\_{\text{\%}}}{f\_{\text{max}}} \tag{12}$$

where *f<sup>e</sup>* and *fmax* are the estimated force and the maximum force of muscle, calculated by Equations (3) and (9), respectively. The figure also shows the three characteristic points defined in Figure 2 to compare the muscle activities with the movements of the hyoid bone and thyroid cartilage. Figure 7 shows the onset and duration of the muscle activity, which was derived from Figure 6. The activation level was defined as values higher than 1% of the maximum force, and the contrast coloring changed depending on the muscle force. The results suggest which muscle mainly contributes to swallowing, and they confirmed that the activities of the muscles correlate with the movements of the hyoid bone and thyroid cartilage.

Figures 2, 6 and 7 show that all suprahyoid muscles (PMH, AMH, GH, AD, PD, and SH) and an upper part of the middle pharyngeal constrictor muscle (PG1) begin to generate contractile forces to pull the hyoid bone in the postero-superior direction in the first phase between *P*<sup>s</sup> and *P*1. Here, ST also contributes a counterbalancing force to the thyroid cartilage, preventing it from being pulled in the posterior direction by the passive force of TH and thyrohyoid membrane, which are stretched. SP seems to contribute to maintaining the posture of thyroid cartilage in this period. In the second phase, between *P*<sup>1</sup> and *P*2, GH, AD, PD, and SH generate larger forces to move the hyoid in an antero-superior direction. Here, TH, SP, and a part of inferior pharyngeal constrictor (PT1 and PT2) begin to generate force to pull the thyroid upward and slightly forward. After reaching the maximum displacement at 1.73 s, in the last phase between *P*<sup>2</sup> and *P*e, the later pharyngeal phase, the activation of all muscles are declined, and the hyoid bone and thyroid cartilage return to the resting position by the passive forces of the stretching muscles and membranes.

௫

%௫

%௫ =

 ௫

= 0

 ௫

**Figure 6.** Estimated results of time changes in muscle forces.

1% 15% 30%~

**Figure 7.** Onsets and durations of activity of swallowing-related muscles.

#### *3.2. Force Direction of the Muscles*

The hyoid bone and thyroid cartilage mainly move in anterior and superior directions, as shown in Figure 2. Therefore, to calculate the maximum contraction forces of the muscles in the anterior-posterior direction and the superior-inferior direction, it is important to investigate the role of each muscle in the movements of the hyoid bone and thyroid cartilage. Figure 8 shows the direction of the maximum contractile force of each muscle in the anterior-posterior and superior-inferior directions, which were calculated from Figure 6.

From Figure 8, the contribution is large in order of AD and GH in pulling the hyoid bone in the anterior direction. For thyroid cartilage, only TH contributes in the anterior movement. In the same manner, the contribution is large in order of PD, SH, PMH, AMH, GH, and AD in elevating the hyoid bone. For the superior movement of thyroid, the contribution order is TH, SP, and PT1.

**Figure 8.** Ratios and directions of maximum force of each muscle. (**a**) Anterior-posterior, (**b**) superior-inferior.

#### **4. Discussion**

#### *4.1. Activity Pattern of the Muscles*

The major contribution of this study was to estimate the sequential activation pattern of the many muscles involved in swallowing based on the dynamics of swallowing, which is hard with conventional methods. In the same manner, the estimated activation pattern can be divided into three phases, corresponding to the sequence of events in swallowing. In the first time period, *P*s–*P*1, the later oral phase, most of the suprahyoid muscles and a part of middle pharyngeal constrictor muscles are activated to elevate the hyoid bone to push the tongue against the hard palate to hold bolus. The result indicates that PMH, AMH, GH, AD, PD, and SH contribute to elevating the hyoid bone and, particularly, PD and SH make significant effort. Subsequently, in the second period, *P*1–*P*2, the early pharyngeal phase, those muscles are more activated to move the hyoid bone anteriorly and superiorly, in which GH and AD show a remarkable contribution to the anterior movement. In addition, TH, SP, and PT1 are activated to raise thyroid cartilage slightly behind the hyoid bone. In the third period, *P*2–*P*e, the later pharyngeal phase, all the muscles tend to relax, returning the hyoid bone and thyroid cartilage to resting positions. During the pharyngeal phase, the larynx rises along with the antero-superior movements of the hyoid bone and thyroid cartilage, and the epiglottis folds backward to cover the entrance of the larynx to keep bolus from entering windpipe. At the same time, the pharynx would extend to transport bolus into the esophagus.

To estimate the muscle activity, electromyography (EMG) was conventionally used as a reliable technique. In particular, fine wire EMG was widely used to evaluate swallowing-related muscles individually and to investigate sequential muscle activation patterns [43–48]. In these previous studies, activities of typical muscles involved in swallowing were investigated, but not the deep muscles, which are difficult to approach by wire electrode. The following findings were obtained:


Similar findings were also observed in this study. That is, GH and AD are activated at an early stage in addition to MH (AMH, PMH), PD, and SH to move the hyoid bone anteriorly and superiorly. The onset of MH (AMH, PMH) precedes the activation of TH, which mainly contributes to the movement of thyroid cartilage, and the activity of SN is quite low. Besides, CP (CP1, CP2) and the lower region of PT (PT3, PT4) slightly exhibited twice activation. However, note that the number of muscles that can be measured by conventional methods are limited, unlike our approach.

#### *4.2. Role of the Muscles*

Figure 8 reveals the role of each muscle in the movements of the hyoid bone and thyroid cartilage. The figure shows the reasonable result, that is, all suprahyoid muscles (AMH, PMH, GH, AD, PD, and SH) contribute to the superior movement of the hyoid bone because they are located above the hyoid bone. Similarly, it is understandable that GH and AD show the most effort in the anterior movement of the hyoid bone as they are in front of the hyoid bone. TH dominantly works on the thyroid movements and it is estimated to exert the largest force in the anterior-superior movement of thyroid cartilage. Meanwhile, TH also contributes to depressing the hyoid bone as it is located between the hyoid bone and thyroid cartilage.

To understand how the muscle morphological properties influence hyoid movement, Pearson et al. [41] investigated muscle attachment sites, and physiological cross-sectional areas (PCSA) of the suprahyoid muscle subsamples (the geniohyoid muscle, stylohyoid muscle, digastric muscle, and mylohyoid muscle) were investigated by an anatomical technique. They evaluated the potential effect of individual muscles on the hyoid bone using the PCSA forces of the muscles. The PCSA force, which represents the strength of the muscles in the anterior-posterior direction or the superior-inferior direction, was calculated from the attachment sites and PCSA of the muscles. As the result, it predicted that GH and MH have the most potential to affect the anterior and superior displacement of the hyoid bone, respectively, which is partially consistent with the findings derived from this study. In addition, our research emphasizes the importance of suprahyoid muscles in terms of swallowing dynamics.

#### *4.3. Limitations*

The findings of this study considering the dynamics of swallowing were consistent with the results of previous studies based on physiological measurement and structural analysis. However, there were some limitations to this study. First, the analysis was of only one healthy young volunteer. Even though the movement of the hyoid bone of the volunteer is a typical triangular trajectory, to generalize the results it would be better to increase the number of volunteers of various ages and genders because hyoid movement is slightly varied among individuals. In particular, to investigate the effects of aging is important, as muscle strength weakens with age and swallowing disorders are common in the elderly. Second, in relation to the first, we developed the musculoskeletal model using the CT data of the young volunteer, but the skeleton geometry also differs among individuals and that difference influences the dynamics calculation. Thus, a deformable standard skeleton model is required to represent different genders and ages. Besides, to analyze the sensitivity of subject-specific modeling to the uncertainties in the identification of muscle attachment site, muscle-tendon parameters and stiffness parameters of soft tissues are interesting because the geometry of head and neck parts is varied among individuals. Third, there are some simplifications. The effect of tongue deformation done by tongue muscles was not accurately simulated but it actually affects the kinematics and posture of the hyoid bone. Similarly, linear elastic and viscous behaviors were assumed for soft tissues surrounding the hyoid bone and thyroid cartilage. Thus, nonhomogeneity and anisotropy of their restrictions were not taken into account. In addition, the inertial effect of the muscles on the movements of the hyoid and thyroid was not considered, because Hill-type muscle model does not include muscle mass. However, the inertia of the attached muscles cannot be ignored as their masses are actually larger than that of the hyoid bone and thyroid. To overcome these simplifications, a 3D, finite element model of muscle and soft tissue should be constructed based on physiological literature. The 3D musculoskeletal model would also enable us to take into account the interaction between muscles and the effect of gravity on muscle deformation. Fourth, there are some assumptions and hypothesis due to the lack of physiological data. Most of the PCSA of the muscles derived from cadaver specimens was used as that of the young volunteer and the PCSA of the pharyngeal constrictor muscles were estimated from a 3D anatomical model. In addition, the estimated muscle forces in this study may not be a unique combination, because the estimated result depended on the objective function of the optimization process. Although the objective function considered physiological properties of muscles and is widely accepted in musculoskeletal studies, its general applicability is not fully accepted. Thus, the results of this study should be validated further by physiological manner.

To investigate muscle activities in swallowing is difficult because swallowing is a complex behavior involving volitional and reflexive activities of many nerves and muscles. In general, EMG measurements are recognized as a reliable method to estimate muscle activities. At the same time, however, there are limitations in measuring deep muscles in the neck even with fine wire electrodes. Further, EMG signals do not necessarily reflect muscle strength. In comparing with EMG measurement, the advantage of the proposed method is that it is able to estimate muscle strengths that can satisfy dynamic consistency parameters with the movement of the hyoid bone and thyroid cartilage calculated using VF images, and as the activities of deep muscles can be estimated. As a result, this technique is only one method to investigate individual muscles in swallowing. If the musculoskeletal analysis of swallowing were validated physiologically, this technique would be useful for swallowing rehabilitation. For example, it may extract the target muscle group, which has to be recovered to obtain the desired therapeutic effect. Another possible application is to propose effective postural compensation maneuver for swallowing safely without aspiration.

#### **5. Conclusions**

A musculoskeletal model was developed, with the main muscles around the hyoid bone and thyroid cartilage included, to study the muscle activity during swallowing. The forces of individual muscles between the later oral phase and the pharyngeal phase of normal swallowing of a healthy volunteer were estimated by inverse dynamics computations with a static optimization technique. The predicted muscle activation patterns showed qualitative agreement with salient features obtained from previous studies. This study shows that the musculoskeletal analysis has the potential to provide deeper insight into swallowing mechanism than conventional methods in terms of swallowing dynamics.

One of our future works is to implement a muscle model with volume like finite element method (FEM) to investigate the effects of deformation of muscle, tongue, and soft tissues on the movements of the hyoid and thyroid. Although the Hill-type muscle model does not have volume, the effect of muscle deformation cannot be ignored in swallowing as the hyoid and thyroid are held by the muscles. Besides, to integrate simulation models of tongue deformation [29,30] and mastication [49], it is also important to analyze the coordination of mastication and swallowing.

**Author Contributions:** Conceptualization, T.H., T.K. (Takahiro Kikuchi), and Y.M.; methodology, T.H. and A.M.; software, T.H., A.M., F.C.-S., and M.U.; data acquisition, T.K. (Takahiro Kikuchi) and Y.M.; writing—original draft preparation, T.H.; writing—review and editing, M.U., T.K. (Takahiro Kikuchi), and Y.M.; supervision, T.K. (Takuji Koike); funding acquisition, T.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by JSPS KAKENHI Grant Numbers 18K09706 and 19K22980 and Tateishi Science and Technology Foundation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Assessment of Motion and Muscle Activation Impacts on Low Back Pain during Pregnancy Using an Inertial Measurement Unit**

**Saori Morino 1,2,\* , Mamoru Yamashita <sup>3</sup> , Fumiko Umezaki <sup>3</sup> , Hiroko Hatanaka <sup>3</sup> and Masaki Takahashi <sup>4</sup>**


Received: 1 May 2020; Accepted: 25 May 2020; Published: 26 May 2020

**Abstract:** Specific physiological changes during pregnancy exert excessive strain on muscles such as the erector spinae (ES) and contribute to low back pain (LBP). The link between LBP and sit-to-stand (STS) motion has previously been investigated through motion analysis using an inertial measurement unit (IMU); however, the factors leading to LBP have not been revealed. Moreover, clinicians require an effective assessment method for reducing the physical burden on pregnant women. Therefore, the investigation of the relationships between motion, muscle load calculated from musculoskeletal model for pregnancy, and the severity of LBP during STS in pregnant women was conducted. Furthermore, this study proposes a method for assessing motion and muscle load during STS using an IMU. The relationship among (i) motion evaluation indices and ES muscle torque, and (ii) the ES torque and the intensity of LBP during STS was investigated. As the results, significant positive correlations were observed between (i) the angular velocity of the torso in the sagittal plane and ES torque, and (ii) two types of evaluation indices of ES torque and intensity of LBP. The proposed method by an IMU attached to the torso could effectively assess ES load related to LBP during STS in pregnant women.

**Keywords:** musculoskeletal model; pregnancy; motion analysis; muscle load evaluation; low back pain

#### **1. Introduction**

Lumbopelvic pain (LPP), which includes low back pain (LBP) and pelvic girdle pain, are common discomforts experienced by women during pregnancy [1]. Symptoms can have a negative impact on the quality of life of the mother in both prenatal and postnatal stages [2,3]. Body weight gain, especially in the abdominal area, and a shift in the position of the body's center of gravity due to fetal growth are distinct changes experienced during pregnancy and a common risk factor of LPP. In addition, an expanded gravid uterus can stretch and weaken the abdominal muscles [4]. These changes in the physical and musculoskeletal system affect the posture and movement of pregnant women and impose strain on muscles and joints, which contributes to pain at various locations such as the lower back [5,6]. For example, the sit-to-stand (STS) motion is a cause of LBP during pregnancy [7]. Hence, physical

methods are required to manage musculoskeletal disturbances during pregnancy, together with the coaching of proper movement to avoid physical loading [8].

Clinical motion analysis typically employs visual observation as the dominant means of assessment; however, more quantitative and objective motion assessments are required to obtain meaningful results [9]. For example, the inertial measurement unit (IMU) allows objective motion measurement without limitations related to the measurement environment or disturbances in motion. Therefore, the IMU has been widely used to analyze motions including gait in a straight path and STS [10,11]. Our previous research investigated the differences in torso motion between pregnant women with and without LPP during STS using an IMU [12]. In the study, an IMU was attached to the lower torso, and angular velocity data in the sagittal plane were used to assess the flexion and extension movements of the trunk. Some key indices were then proposed, that is, maximum peak, minimum peak, and peak-to-peak (PP). The maximum peak represents forward motion of the torso, and the minimum peak represents backward motion, with larger values indicating faster motion. PP is the difference between the maximum and minimum peaks and therefore represents a shift in motion from forward to backward. The results of the study revealed that the maximum peak was greater in the LPP group than in the non-LPP group. However, the role of other motion characteristics as risk factors of pain, such as additional stress on lumbar muscles due to muscle load, remains unclear.

Physical movement therapy based on assessments of motion analysis and muscle load, such as the exerted muscle force, increases the effectiveness of physical therapy for relieving pain. Evaluation of the erector spinae (ES) muscle force is particularly useful for managing LBP because muscle load and fatigue are two of the main risk factors associated with LBP during pregnancy [1]. Although the musculoskeletal model can be used to estimate muscle force, it is difficult to estimate the co-contraction activation of agonist and antagonist muscles, such as the ES and rectus abdominis (RA) [13,14]. To address these issues, we previously proposed a method for estimating muscle torque from the musculoskeletal model using a genetic algorithm (GA) for pregnant women [15].

STS motion assessment should consider the ES muscle load in order to manage LBP in pregnant women. However, the link between muscle load during STS motion and LBP has not yet been revealed for pregnant women. Therefore, the aim of this study was to propose a method for assessing motion and muscle load during STS motion using an IMU. A secondary aim was to investigate the relationship among STS motion, muscle load, and LBP during pregnancy using the proposed assessment method.

#### **2. Materials and Methods**

#### *2.1. Study Design*

This is a cross-sectional study. Body motion analysis of STS motion in pregnant women was performed to investigate the relationship among STS motion, muscle load, and LBP. The motion condition was determined using pitch angular velocity data obtained from an IMU attached to the lower trunk of pregnant participants with LBP to evaluate their specific motions [16]. Subsequently, musculoskeletal models tailored to each participant in the study were constructed, and the muscle torque was estimated using a GA. Finally, the different relationships were investigated using the value of motion indices, estimated ES muscle torque, and intensity of LBP during pregnancy. The STROBE checklist of this study is Table S1.

#### *2.2. Participants*

The study protocol was approved by the Ethics Committee (Kishokai Medical Corporation: approval no. 2015\_002). Pregnant women who were undergoing a prenatal health checkup in an obstetrics and gynecology clinic were invited to participate in this study. The purpose and protocol of this study were explained to all participants by the examiner. Written informed consent was obtained from each participant, including consent to participate and to publish the findings. Few studies have focused on the biomechanics of pregnant women and, moreover, studies that investigated the

relationships between biomechanical factor and some complaints by motion analysis for pregnancy could not be found. Then, some studies that focused on the biomechanics of pregnant women in Japan were referenced for determining sample size [17,18]. In each study, eight pregnant women were investigated. Therefore, the sample size of this study was decided as the same number of these previous studies, and thus 11 women were recruited, taking account of the case of dropout. Hence, a group of 11 pregnant women participated in this study. The average duration of pregnancy was 34.5 weeks (standard deviation = 2.0 weeks), and all participants had singleton pregnancies. None of the participants suffered from any serious neurological or orthopedic conditions, particularly those associated with LBP, such as intervertebral disk displacement or spinal cord injury. Before the STS motion analysis, it was confirmed that the participants had not suffered any external injuries that could influence the analysis. Additionally, physical characteristics (age, height, body mass at the time of experiment, body mass gained during pregnancy) were obtained via a questionnaire. Personal data of the participants are listed in Table 1.



<sup>1</sup> LBP: low back pain; Values are mean values (standard deviation) with [range], except for the number of women with LBP before pregnancy. *p*-values are the result of the Mann–Whitney *U* test.

#### *2.3. Sensors for the STS Motion Experiment*

Motion, force, and electromyographic (EMG) signal data were obtained using the protocol of our previous study to calculate joint torque and estimate muscle torque [15]. An IMU (TSND151, ATR-Promotions Co., Ltd., Kyoto, Japan) was used to obtain joint angles and segment orientations during the STS motion analysis. The IMU contains tri-axial gyroscopes, accelerometers, and magnetometers and has been used for human motion analysis [19,20]. In addition, the IMU can be used in synchronization with the EMG electrodes that were used in this study. The participants were assessed using six IMUs. Five IMUs were bilaterally attached to a fixed belt, 10 cm above the lateral malleolus, 10 cm above the patella, and at the level of the L4 spinous process in order to measure the bilateral motion of the shank segments, thigh segment, and lower trunk (pelvis), respectively (Figure 1a). Another IMU was attached to the level of the C7 spinous process using skin tape to measure the motion of the upper trunk. A Nintendo Wii balance board (WBB; Nintendo, Kyoto, Japan) was used to measure the vertical reaction forces (Figure 1b). WBB has been recently proposed as an easily available device for measuring ground reaction force and center of pressure (CoP) displacement and an inexpensive alternative to force plates, which are common in laboratories for the assessment of posture and balance [21,22]. To measure the reaction force from the chair and the ground reaction force, respectively, one WBB was placed on the chair, and the other was placed under the participant's feet. Four surface EMG electrodes (SE-C-AMP-H100, ATR-Promotions Co., Ltd., Kyoto, Japan) were used to measure muscle activation of the ES and RA muscles (Figure 1a). Electromyography electrodes

were bilaterally placed in parallel with the fiber orientation of the RA and longissimus (part of the ES) muscles.

**Figure 1.** Experimental settings: (**a**) position of the inertial measurement units (IMUs) and surface electrodes; (**b**) experimental settings.

#### *2.4. Participant-Specific Musculoskeletal Models*

Full-body musculoskeletal models were constructed for each pregnant woman using Biomechanics of Bodies (BoB; Marlbrook Ltd., Bromsgrove, United Kingdom). BoB is a package containing the musculoskeletal model of a human being, 36 skeletal segments, and over 600 muscle units [23]. BoB's biomechanical software is based on the MATLAB and Simulink environments [24]. Torque in the joints corresponding to the observed motion and force data was calculated by inverse dynamic analysis using BoB, and the muscle torque was estimated using an optimization approach [25,26]. The following process was then undertaken to develop musculoskeletal models for each subject in accordance with our previously proposed procedure [15].

A skeletal model previously modified to the average body type of Japanese women was used [15]. Then, the total body height and weight were adjusted for each participant. Furthermore, the muscular model was also changed because the abdominal muscles were stretched and weakened during pregnancy by an enlarged gravid uterus [15]. In the previous study, a model of the stretched RA muscle was constructed using information from an anthropometric dimensions database for pregnant women and measurements of the three height levels of waist girth for the participants. Thus, the same three height levels of waist girth were measured for each participant in this study. Weight gain, which is a characteristic change during pregnancy, was also considered for each participant.

#### *2.5. Experimental Settings and STS Motion Conditions*

An overview of the experimental settings is given in Figure 1b. For the motion task, participants with IMUs attached to their body segments were asked to start from a standing position in front of a chair. Then, they were instructed to sit down and stand up using the chair as a support. On the basis of the results of our previous study [12], two types of STS motion, characterized by different values of maximum peak angular velocities of the torso in the sagittal plane (slow and fast), were conducted three times by each participant. The difference between the two types of STS motion was the forward incline speed of participant trunks. The incline speed was checked and instructed by the examiner. The participants moved at the prescribed velocities and performed the STS motion without using hands. The participants were asked to put their arms along their sides. The maximum peak represents the maximum velocity of forward trunk movement during STS motion (an example of the indices from

our previous research is shown in Figure A1). Specifically, the participants were instructed to move at low angular velocities for "trial slow" and at high angular velocities for "trial fast", so that the IMU attached to the lower trunk would indicate a value similar to that given in Table A1. The order of the type of motion was determined randomly. A rest period of approximately 10 s was fixed between each motion. An armless, backless chair was used and adjusted to the height of the participant's knee in the condition with WBB, in accordance with previous research [27]. The STS movement was performed by the participants without their shoes on, and the medial borders of their feet were placed 10–15 cm apart.

#### *2.6. Assessment of LBP*

During the first visit, the occurrence of LBP in the year prior to the pregnancy was investigated. Then, the participants were asked about the presence or absence of LBP during the STS motion analysis in every trial. If the participants reported LBP, the pain intensity was assessed using an 11-point numerical rating scale (NRS) [28], where the low and high endpoints represent the extremes of "no pain" and "worst pain", respectively. On the basis of the ratings, NRS > 0 was defined as the condition for the presence of LBP. The NRS score during STS motion was obtained for each participant.

#### *2.7. Processing of Measured Data*

Signal processing of IMU data during STS motion analysis was performed using MATLAB Release 2016a (MathWorks, Release 2016a, Natick, MA, USA). The coordinate system in the experimental settings is shown in Figure 2. First, the segment orientation and joint angles of each segment were calculated as motion data using IMU data and the same protocol as our previous study [15]. The calculated orientations of body segments and joint angles are shown in Figure 2. The direction of the axis of each IMU was determined as the line of each segment.

**Figure 2.** Coordinate system (**left**) and definition of orientations and angles (**right**). Calculated orientations of body segments and joint angles from each IMU are shown, where θ and φ correspond to the symbols in Equation (1).

The joint angle was defined as 0◦ for an anatomically upright posture. The vertical axis in the sagittal plane was used as the reference axis; thus, it was defined as 0◦ to calculate the inclination angle. Joint angles were calculated using Equation (1), where θ and φ correspond to the symbols used in Figure 2.

$$\begin{aligned} \theta\_1 &= \varphi\_1 \\ \theta\_2 &= 180 + \varphi\_2 - \varphi\_1 \\ \theta\_3 &= \varphi\_3 - \varphi\_2 \\ \theta\_4 &= \varphi\_4 - \varphi\_3 \end{aligned} \tag{1}$$

*τ τ*

௦௧ሺሻ = ଵሺሻ௫<sup>మ</sup> + ଷሺሻ௫<sup>ర</sup> *τ* Next, the ground reaction forces were calculated using data from the WBBs (body mass acting on the WBBs, kg) under the foot and hip. Thus, the force (N) in each trial was calculated using Equation (2). Then, these motion and force data were inputted to the BoB package with the modified musculoskeletal

௦௧\_ாௌሺሻ = ହாௌሺሻ௫<sup>ల</sup> + ாௌሺሻ௫<sup>ఴ</sup>

௦௧ሺሻ = ௦௧\_ோሺሻ + ௦௧\_ாௌሺሻ

1 ∑ ሺ௧ሺሻ − ௦௧ሺሻሻ <sup>ଶ</sup>

*τ τ*

ሺ௧ሺሻ − ௦௧ሺሻሻ<sup>ଶ</sup>

ிௌ =

ோெௌா = ඩ

*τ*

ୀଵ

1  ୀଵ

models for pregnant women in order to calculate the joint torque. Finally, the optimization process in BoB was performed using the joint torque to estimate the muscle torque in the model.

$$F = ma\tag{2}$$

Here, *F* is the reaction force (N), *m* is the body mass measured by the WBB (kg), and *a* is the acceleration due to gravity (m/s 2 ).

#### *2.8. Co-Contraction Activation Estimation Using GA*

Muscle torque as an index of muscle load was calculated from obtained data in this study. An overview of the system used to estimate muscle torque is presented in Figure 3. For each trial, muscle torque was estimated by considering co-contraction activation using a GA according to the following process. The procedure for estimating RA and ES muscle torque is fundamentally the same as that in our previous study [15]. First, Equation (3) was used to estimate the muscle torque. The equation is a mathematical model proposed by Oyong and Jauw [29] that converts EMG signals to muscle torque:

$$
\pi\_{\rm est}(k) = \mathfrak{x}\_1 E(k)^{\mathfrak{x}\_2} + \mathfrak{x}\_3 E(k)^{\mathfrak{x}\_4} \tag{3}
$$

where *k* is the sampling time, τ*est* (*k*) is the estimated torque, *E*(*k*) is the processed EMG signal, and *x<sup>j</sup>* represents the associated model parameters. Accordingly, the torque of each RA and ES muscle, as well as the sum of the two muscle torque results, were calculated using Equation (4): *τ τ*

$$\begin{aligned} \tau\_{\text{est\\_}}(k) &= \mathbf{x}\_1 E\_{RA}(k)^{\mathbf{x}\_2} + \mathbf{x}\_3 E\_{RA}(k)^{\mathbf{x}\_4} \\ \tau\_{\text{est\\_}ES}(k) &= \mathbf{x}\_5 E\_{ES}(k)^{\mathbf{x}\_6} + \mathbf{x}\_7 E\_{ES}(k)^{\mathbf{x}\_8} \\ \tau\_{\text{est}}(k) &= \tau\_{\text{est\\_}RA}(k) + \tau\_{\text{est\\_}ES}(k) \end{aligned} \tag{4}$$

where *ERA* is the processed EMG signal of RA, and *EES* is the processed EMG signal of ES. *τ τ*

**Figure 3.** Overview of the system used to estimate muscle torque considering co-contraction activation using a genetic algorithm (GA). EMG: electromyography, IMU: inertial measurement unit.

Trial conditions in the GA process such as population size were the same as those in our previous study [15]. Subsequently, the joint torque during STS motion was calculated using the motion, force, and skeletal data via inverse dynamic calculation in the pregnant musculoskeletal model, and the

ாௌ\_ோெௌ = ඩ

1  ୀଵ

ሺ௦௧\_ாௌሺሻሻ<sup>ଶ</sup>

value of torque was used as the actual torque in the GA. These values enabled us to estimate the muscle torque via an optimized GA using Equation (5). The fitness values (*VFS*) in this equation can be maximized only when a certain estimated muscle torque (τ*est*) fits the actual joint torque (τ*act*) as a supervisory signal.

$$V\_{FS} = \frac{1}{\sum\_{k=1}^{n} \left(\tau\_{act}(k) - \tau\_{est}(k)\right)^2} \tag{5}$$

Here, *VFS* is the fitness value, *n* is the number of τ*act* and τ*est* values in each experiment trial, and τ*act* is the joint torque calculated from the musculoskeletal model. In addition, a completion condition was set in this study. To assess the precision of the estimated muscle torque, the root mean square error (RMSE) was calculated using Equation (6):

$$d\_{RMSE} = \sqrt{\frac{1}{n} \sum\_{k=1}^{n} \left(\tau\_{act}(k) - \tau\_{est}(k)\right)^2} \tag{6}$$

where τ*act* is the actual joint torque during movement, and τ*est* is the estimated muscle torque. Smaller RMSE values indicate that the muscle torque values were better estimated. From a previous study that estimated lumbar muscle force using the EMG-based model [30], the mean ratio of RMSE between the measured and estimated lumbar sagittal moments to the maximum peak of the measured sagittal moments was approximately 22.34% in the sagittal plane. In the study, RMSE values that were higher than those obtained by previous research were not considered suitable for estimated muscle torque [31,32]. Thus, the cut-off completion condition value in this study was 20% of the peak value of τ*act* in each trial. The GA step was completed when the value of RMSE was less than 20% of the τ*act* value in each trial.

#### *2.9. Evaluation Index of Muscle Torque*

A universal evaluation index calculated from muscle torque during STS motion has not been established by previous research. Thus, the following indices were calculated to investigate the relationship between the motion features and muscle torque of the ES during STS motion in pregnant women (Figure 4). First, the value of the ES muscle torque (MT) during STS motion at the maximum peak of the angular velocity in the sagittal plane of the lower trunk segment (MP) was detected and termed MT-MP (Figure 4, point A). Then, the difference in muscle torque at MP and the standing position (DMT-MP-stand) was also calculated (Figure 4, point B) to remove the individual differences in muscle torque in the standing position. The maximum peak of the ES muscle torque during STS motion (Max-MT) was also detected (Figure 4, point C), and the difference between the value of ES muscle torque during STS motion and the standing position (DMT-Max-stand) was calculated (Figure 4, point D). Furthermore, the root mean square (RMS-MT), which is commonly used to express the effective value of a waveform, was calculated using the ES muscle torque data for each motion by Equation (7). All relevant indices were calculated for each trial and each participant.

$$d\_{ES\\_RMS} = \sqrt{\frac{1}{n} \sum\_{k=1}^{n} \left(\tau\_{est\\_ES}(k)\right)^2} \tag{7}$$

**Figure 4.** Indices for muscle torque evaluation.

#### *2.10. Statistical Analysis to Compare Muscle Torque Evaluation Indices*

The relationship among the indices of angular velocity, which represents motion information, the estimated ES muscle torque, and LBP, was investigated using participant data from the motion analysis. Before the investigation, a normality of each of the data was checked by using of a Shapiro–Wilk test. Then, a non-parametric test was conducted because none of the datasets followed a normal distribution. First, a Wilcoxon signed-rank test was conducted to investigate the differences among the five indices of the estimated ES muscle torque (MT-MP, DMT-MP-stand, Max-MT, DMT-Max-stand, and RMS-MT) between the STS motion in "trial slow" and "trial fast". That is, the five indices were compared for all 11 participants in the two trials.

Next, a bivariate correlation analysis using the Spearman's rank correlation coefficient was conducted to investigate the correlation between the five indices of estimated ES muscle torque and the maximum peak during STS motion. In this analysis, PP was investigated instead of the maximum peak because the timing of Max-MT was almost the same as PP among the motion indices, as observed in Figures 4 and A1. Moreover, the data from "Trial fast" representing the characteristic motion of pregnant women who were experiencing pain were used to investigate the motions that could be linked to LBP. Furthermore, Spearman's rank correlation coefficients between each index of the estimated ES muscle torque and the maximum peak during STS motion, as well as between each index and PP value during STS motion, were calculated. Finally, a Spearman's rank correlation coefficient was calculated using five indices of the estimated ES muscle torque and the LBP-related NRS score to investigate the correlation between the degree of muscle load and the intensity of LBP. Nine pregnant women reported LBP during "trial fast," but none reported LBP during "trial slow". Thus, data of these nine participants were used in the last statistical analysis. The analyses were performed on SPSS 23.0 (IBM SPSS Statistics, Chicago, IL, USA) with a significance threshold of *p* > 0.05.

#### **3. Results**

#### *3.1. Di*ff*erence of Muscle Torque According to Di*ff*erence of Motion Trial*

The results of the Wilcoxon signed-rank test are presented in Table 2. The DMT-Max-stand index of "trial fast" was significantly higher than that of "trial slow." No significant differences were observed for the other muscle torque indexes between each trial.

**Table 2.** Differences in evaluated indices of muscle torque according to type of sit-to-stand (STS) motion (*N* = 11).


<sup>1</sup> MT-MP: the value of the erector spinae (ES) muscle torque (MT) during STS motion at the maximum peak (MP) of the angular velocity; <sup>2</sup> DMT-MP-stand: the difference in muscle torque at MP and the standing position; <sup>3</sup> Max-MT: the maximum peak of the ES muscle torque during STS motion; <sup>4</sup> DMT-Max-stand: the difference between the value of ES muscle torque during STS motion and the standing position; <sup>5</sup> RMS-MT: the root mean square of the ES muscle torque. Values shown are mean values (standard deviation). The data in bold is statistically significant.

#### *3.2. Correlation between Indices of Estimated Muscle Torque and the Value of Motion Indices in STS*

Table 3 lists the results of the statistical correlation analysis performed for the five indices of the estimated ES muscle torque and values of the motion index, revealing that a significant positive correlation was observed between DMT-Max-stand and PP. No significant correlations were observed between maximum peak and indices of estimated muscle torque.


RMS-MT <sup>6</sup> 0.112 0.536

Max-MT 0.282 0.112 **DMT-Max-stand 0.355 0.043** MT-MP 0.141 0.434 DMT-MP-stand 0.145 0.419 RMS-MT 0.176 0.327

**Table 3.** Correlation analysis results for motion and muscle torque (*N* = 11).

1 *r*: correlation coefficient; <sup>2</sup> Max-MT: the maximum peak of the ES muscle torque during STS motion; <sup>3</sup> DMT-Max-stand: the difference between the value of ES muscle torque during STS motion and the standing position; <sup>4</sup> MT-MP: the value of the ES muscle torque (MT) during STS motion at the maximum peak (MP) of the angular velocity; <sup>5</sup> DMT-MP-stand: the difference in muscle torque at MP and the standing position; <sup>6</sup> RMS-MT: the root mean square of the ES muscle torque; <sup>7</sup> PP: peak-to-peak. The data in bold is statistically significant.

#### *3.3. Correlation between the Indices of Estimated Muscle Torque and the Intensity of LBP in STS*

PP <sup>7</sup>

Table 4 presents the results of the statistical correlation analysis for the five indices of estimated ES muscle torque and the LBP-related NRS score, revealing significant positive correlations between (i) Max-MT and LBP intensity and (ii) DMT-MP-stand and LBP intensity. No other significant correlations were observed between the indices of estimated muscle torque and LBP intensity.


**Table 4.** Correlation analysis results for muscle torque and low back pain (*n* = 9).

1 *r*: correlation coefficient; <sup>2</sup> Max-MT: the maximum peak of the ES muscle torque during STS motion; <sup>3</sup> DMT-Max-stand: the difference between the value of ES muscle torque during STS motion and the standing position; <sup>4</sup> MT-MP: the value of the ES muscle torque (MT) during STS motion at the maximum peak (MP) of the angular velocity; <sup>5</sup> DMT-MP-stand: the difference in muscle torque at MP and the standing position; <sup>6</sup> RMS-MT: the root mean square of the ES muscle torque. The data in bold is statistically significant.

#### **4. Discussion**

This study proposed a method for assessing motion and muscle load in pregnant women during STS motion on the basis of an inertial measurement unit. This study also investigated the relationships between motion evaluation indices, muscle load (erector spinae activation), and the severity of LBP. The link between muscle load during STS motion and LBP in pregnant women has not previously been determined, despite this being a common discomfort experienced by women during pregnancy, which can have a negative impact on quality of life during both prenatal and postnatal stages. As its results, this study established the relationship between erector spinae muscle torque and LBP severity during STS motion in pregnant women.

#### *4.1. Di*ff*erence of Muscle Torque According to Di*ff*erence of Motion Trial*

The results in Table 2 suggest that the ES muscle torque during STS motion, excluding muscle torque in a standing still position, was higher for specific motions exhibited by pregnant participants with LBP than for those exhibited by participants without LBP. These results indicate that higher muscle load from an immobile standing position could be related to LBP during STS motion.

#### *4.2. Correlation between Indices of Estimated Muscle Torque and the Value of Motion Indices in STS*

The positive correlation between DMT-Max-stand and PP, which is a parameter describing the change in incline speed from forward (maximum peak) to backward (minimum peak) motion, suggests that ES muscle load might increase with a higher degree of shift change from forward bending to backward bending motion. Therefore, the results confirm that the value of maximum ES muscle torque in relation to the motion of standing up can be evaluated using the PP of angular velocity from the IMU attached to the lower torso, without the need for a direct evaluation of muscle load.

#### *4.3. Correlation between the Indices of Estimated Muscle Torque and the Intensity of LBP in STS*

In addition to Table 4, scatter diagrams of the NRS score with the Max-MT and DMT-MP-stand are shown in Figure A2. The significant positive correlation between Max-MT and LBP intensity implies that higher values of maximum ES muscle torque during STS motion could indicate a higher correlation with LBP. Moreover, a significant positive correlation between DMT-MP-stand and LBP intensity was also observed. Thus, a larger value of maximum ES muscle torque in relation to the standing up motion might be related to more severe LBP symptoms.

#### *4.4. Relationships among the Motion Indices from an IMU, Estimated Muscle Torque, and the Intensity of LBP in STS*

Our previous research revealed that the STS motion particular to pregnant women with LBP can be observed using an IMU attached to the lower torso. For example, a higher maximum angular velocity peak in the sagittal plane was observed in the LBP group. However, other risk factors, such as muscle fatigue, were not considered in that research. In the current study, however, the relationships among the motion index, muscle torque, and LBP during STS motion were investigated for pregnant women using motion, force, and EMG data, respectively. The results in Section 3.2 reveal a statistically significant positive correlation between PP and ES muscle torque (out of muscle torque in the standing position) during STS motion. Simultaneously, the results in Section 3.3 reveal a statistically significant positive correlation between the maximum muscle torque and the intensity of LBP. These results suggest that the PP of angular velocity in the sagittal plane, obtained from an IMU attached to the lower trunk, may have expressed not only the trunk motion characteristics, but also the ES muscle torque exhibited during STS motion for pregnant women. Mechanical information, such as muscle torque estimated by the proposed model, is more useful for quantifying the physical burden related to LBP. In addition, the results reveal the relationship between higher ES muscle torque and LBP severity during STS motion in pregnant women. Thus, this study indicates the potential for conducting pregnancy motion analysis and muscle torque evaluation, which may be associated with LBP, using the proposed IMU-based method, without the need for EMG measurement.

This is the first study that has focused on the biomechanics of pregnant women and investigated the relationships between ES muscle torque during motion and LBP by motion analysis. However, there are various factors that should be considered when investigating the relationships between biomechanics and LBP during pregnancy. For example, the iliopsoas muscle has a relevant role in STS movement and is related to LBP. Moreover, this muscle might be stretched as well as the rectus abdominis during pregnancy. Similarly, other muscles such as the diaphragm and other abdominal wall muscles might be related to LBP [33,34]. Moreover, other studies have indicated that there is a relationship between myofascial trigger points in the paraspinal muscles and LBP [35]. Therefore, further study would be required to investigate the existence of an extensive and detailed biomechanical factor related to LBP during pregnancy.

#### **5. Conclusions**

Experiments on STS motion with pregnant women indicated the possibility that muscle torque related to LBP can be evaluated using an IMU. Higher ES muscle load from an immobile standing base state might be related to LBP during STS motion in pregnancy. However, measuring muscle load is difficult, especially during STS motion and for pregnant women. Our results revealed that the PP of angular velocity in the sagittal plane measured by an IMU attached to the lower trunk indicates greater ES muscle torque during STS motion. In addition, our experiments revealed the relationship between higher ES muscle torque and LBP severity during STS motion in pregnancy.

The limitations of this study include the small sample size and cross-sectional design. Therefore, the causal relationship between motion characteristics, exerted ES muscle torque, and LBP during STS motion could not be revealed. Moreover, the value of correlation coefficient of the relationship between PP and DMT-Max-stand was not high. Thus, the correlation might not be strong. Similarly, the results of the relationships between muscle torque and intensity of LBP were not enough to decide the higher muscle torque indicating the higher intensity of LBP. Therefore, other factors should be considered for utilizing the results of this study. Hence, further research that includes a larger number of participants, a longitudinal design, and surveying of other risk factors of LBP is required to support the results of this study. Despite these limitations, this study established the effectiveness of an IMU attached to the torso for assessing STS motion in pregnant women. Moreover, the proposed physical assessment method can be employed to manage LBP during pregnancy.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-3417/10/11/3690/s1, Table S1: STROBE Checklist.

**Author Contributions:** Conceptualization, S.M. and M.T.; methodology, S.M., F.U., and M.Y.; software, S.M.; formal analysis, S.M. and M.T.; investigation, S.M. and H.H.; resources, M.Y., F.U., and M.Y.; data curation, S.M., F.U., and H.H.; writing—original draft preparation, S.M.; writing—review and editing, M.T.; visualization, S.M.; supervision, M.Y.; project administration, M.Y. and M.T.; funding acquisition, S.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by JSPS KAKENHI under grant 20K19348.

**Acknowledgments:** We are grateful to the staff of Kishokai for recruiting the participants and their cooperation in obtaining the assessments. We are also grateful to the pregnant women for their cooperation and participation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Wave forms and indices of angular velocity of lower trunk on the sagittal plane.

**Table A1.** Value of angular velocity of trunk motion in the sagittal plane, according to the occurrence of lumbopelvic pain (LPP) in our previous research.


A Student's *t*-test was used; values are shown as mean (standard deviation); LPP: lumbopelvic pain, PP: peak-to-peak. The data in bold is statistically significant.

**Figure A2.** Scatter diagram of numerical rating scale (NRS) score with maximum of muscle torque and difference of maximum of muscle torque with stand.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1817-6