Subject-specific Musculoskeletal Model for Studying Bone Strain During Dynamic Motion

Summary

During landing, lower-body bones experience large mechanical loads and are deformed. It is essential to measure bone deformation to better understand the mechanisms of bone stress injuries associated with impacts. A novel approach integrating subject-specific musculoskeletal modeling and finite element analysis is used to measure tibial strain during dynamic movements.

Abstract

Bone stress injuries are common in sports and military trainings. Repetitive large ground impact forces during training could be the cause. It is essential to determine the effect of high ground impact forces on lower-body bone deformation to better understand the mechanisms of bone stress injuries. Conventional strain gauge measurement has been used to study in vivo tibia deformation. This method is associated with limitations including invasiveness of the procedure, involvement of few human subjects, and limited strain data from small bone surface areas. The current study intends to introduce a novel approach to study tibia bone strain under high impact loading conditions. A subject-specific musculoskeletal model was created to represent a healthy male (19 years, 80 kg, 1,800 mm). A flexible finite element tibia model was created based on a computed tomography (CT) scan of the subject's right tibia. Laboratory motion capture was performed to obtain kinematics and ground reaction forces of drop-landings from different heights (26, 39, 52 cm). Multibody dynamic computer simulations combined with a modal analysis of the flexible tibia were performed to quantify tibia strain during drop-landings. Calculated tibia strain data were in good agreement with previous in vivo studies. It is evident that this non-invasive approach can be applied to study tibia bone strain during high impact activities for a large cohort, which will lead to a better understanding of injury mechanism of tibia stress fractures.

Introduction

Bone stress injuries, such as stress fractures, are severe overuse injuries requiring long periods of recovery and incurring significant medical cost^1,2. Stress fractures are common both in athletic and military populations. Among all sports related injuries, stress fractures account for 10% of the total³. In particular, track athletes face a higher injury rate at 20%⁴. Soldiers also experience a high rate of stress fractures. For instance, a 6% injury rate was reported for the US Army¹ and a 31% injury rate was reported in the Israeli Army⁵. Among all reported stress fractures, tibia stress fracture is the most common one^6,7,8.

Sports and physical trainings with a higher risk of tibia stress fracture are normally associated with high ground impacts (e.g., jumping, landing, and cutting). During locomotion, a ground impact force is applied to the body when the foot contacts the ground. This impact force is dissipated by the musculoskeletal system and footwear. The skeletal system serves as a series of levers allowing muscles to apply forces to absorb the ground impact⁹. When the leg muscles cannot adequately reduce the ground impact, lower-body bones must absorb the residual force. Bone structure will experience deformation during this process. Repetitive absorption of residual impact force may result in microdamages in the bone, which will accumulate and become stress fractures. To date, information related to bone reaction to external ground impact forces is limited. It is important to study how the tibia bone responds to the mechanical load introduced by high impact forces during dynamic motions. Examining tibia bone deformation during high impact activities could lead to a better understanding of the mechanism of tibia stress fracture.

Conventional techniques used to measure bone deformation in vivo rely on instrumented strain gauges^{10,11,12,13,14,15}. Surgical procedures are needed to implant strain gauges on bone surface. Due to the invasive nature, in vivo studies are limited by a small sample of volunteers. In addition, the strain gauge can only monitor a small region of the bone surface. Recently, a non-invasive method utilizing computer simulation to analyze bone strain was introduced^16,17. This methodology allows for the ability to combine musculoskeletal modeling and computational simulations to study bone strain during human movement.

A musculoskeletal model is represented by a skeleton and skeletal muscles. The skeleton consists of bone segments, which are rigid or non-deformable bodies. Skeletal muscles are modeled as controllers using the progressive-integral-derivative (PID) algorithm. The three-term PID control uses errors in estimation to improve the output accuracy¹⁸. In essence, PID controllers representing muscles try to duplicate body movements by developing necessary forces to produce length changes of the muscles over time. The PID controller uses the error in the length/time curve to modify the force for reproducing the movement. This simulation process creates a feasible solution to coordinate all muscles to work together to move the skeleton and produce the body movement.

One or more segments in the skeleton of the musculoskeletal model can be modeled as flexible bodies to allow measurement of deformation. For instance, the tibia bone can be broken down into a finite number of elements, which consists of thousands of elements and nodes. The effect of mechanical loading on the flexible tibia can be examined through finite element (FE) analysis. The FE analysis calculates the loading response of individual elements over time. As the number of bone elements and nodes increase, the computation time of the FE analysis will significantly increase.

To reduce computational cost with accurate evaluation of flexible bodies' deformation, modal FE analysis has been developed and used within the automotive and aerospace industry^19,20. Instead of analyzing individual FE elements' responses to mechanical load in the time domain, this procedure assesses an object's mechanical responses based on different vibrational frequencies in the frequency domain. This method results in a significant reduction in computation time while providing accurate measurement of deformation²⁰. Although modal FE analysis has been widely used to study mechanical fatigue in automotive and aerospace areas, the application of this method has been very limited in human movement science. Al Nazer et al., used a modal FE analysis to examine tibial deformation during human gait and reported encouraging results^16,17. However, their method was greatly affected by only using limited kinematic data from an experiment to drive the computer simulations; There were no real ground impact forces used to assist the simulations. This approach may be reasonable for studying low impact slow motions such as walking, but it is not a feasible solution to study high ground impact movements. Thus, in order to examine lower-body bone reactions during dynamic high impact activities, it is essential to develop an innovative approach to address the limitations associated with the previously reported method. Specifically, a method utilizing accurate experimental kinematic data and real ground impact forces must be developed. Therefore, the goal of this study was to develop a subject-specific musculoskeletal model to perform multibody dynamic simulations with modal FE analysis to examine tibial strain during high impact activities. A dynamic high impact movement represented by drop-landings from different heights was selected to test the method.

Protocol

The experiment was conducted under the Helsinki Declaration. Prior to data collection, the subject reviewed and signed the consent form approved by the University Institutional Review Board before participating in the study.

1. CT Imaging Protocol

1. Take the participant to a facility where a CT scanner is housed. Prior to the CT scan, configure the CT machine with the following parameters: CT slice thickness of 0.625 mm, CT field of view of 15 cm x 15 cm, and auto setting for parameters of peak kilo-voltage (kVp) and milliampere-seconds (mAs) using machine algorithm.
2. Ask the participant to lie on a table that slides into a ring in the CT scanner. Ask the participant to remain very still during the CT scan. Scan each leg separately from the calcaneus through the distal end of the femur.
3. Upon completion of the CT scan, export the CT images in a digital imaging and communications in medicine (DICOM) format. Choose an image size of 512 x 512 pixels (gray scale).
   NOTE: The CT imaging protocol normally lasts less than 1 h. The radiation dose is minimal. It presents no greater risk than the that encountered during normal X-ray medical procedures.

2. Anthropometric Measurement Protocol

1. During the laboratory visit, prior to motion capture, measure the participant's body mass (kg), body height (mm), distance between the anterior-superior iliac spines (ASISs) (mm), leg length (mm), knee joint width (mm), and ankle joint width (mm).
2. Inter ASIS distance measurement: Use a caliper to measure the linear distance between the left ASIS and right ASIS.
3. Leg length measurement: Use a tape measure to measure the linear distance the ASIS and medial malleolus for both legs.
4. Knee joint width measurement: Use a caliper to measure the linear distance between the lateral and medial epicondyles of the femur for both knees.
5. Ankle joint width measurement: Use a caliper to measure the linear distance between the lateral and medial malleoli for both legs.
   NOTE: The inter ASIS distance, leg length, knee and ankle width are used to build a subject model in a biomechanics software (see Table of Materials) for performing kinematic and kinetic calculations.

3. Motion Capture Protocol

NOTE: See Table of Materials for all software and tools used.

1. Placement of reflective markers
   1. Place 14-mm reflective markers on the participant's body at the following anatomical bony landmarks: acromion processes, sternoclavicular joints, base of sternum, posterior process of the 10^th thoracic vertebrae, ASISs, posterior-superior iliac spines (PSISs), 1.5 cm above the lateral knee joint lines, 1.5 cm above the medial knee joint lines, lateral malleoli, medial malleoli, posterior heels, bases of the second metatarsals, and bases of the fifth metatarsals.
   2. Place semi-rigid plastic plates with 4-marker clusters on the thighs and shanks, respectively.
      NOTE: For a better motion capture outcome, the participant is recommended to be barefoot and wear skin-tight clothing. In addition, the marker placement procedure follows a modified "Plug-in-Gait" protocol²¹. A total of 39 reflective markers are used for the motion capture and 34 of them are attached to the lower body.
2. Instruct the participant to warm up by walking on a motorized treadmill at a self-selected speed for 5 min.
3. Calibration of room space for motion capture procedure
   1. Power on the motion capture system (12 high-speed infra-red cameras) and two force plates. Open a motion capture software program. Within the main program window, open the 'Resources' pane. Click the "System" tab. Configure the camera frequency at 200 Hz and force plate frequency at 2,000 Hz.
   2. Within the main program window, open the 'Tools' pane. Click the "System Preparation" button. Click "Calibrate Cameras". Click "Start". Ask a research staff to wave a standard 5-marker calibration wand to perform a dynamic calibration within the room space where the drop-landing movements are to be performed. Click "Stop" after 5 s of wand data have been acquired.
   3. Place the calibration wand flat on the floor to align with a corner of a force plate for the purpose of specifying a reference location (origin) for the calibrated space. Click "Set Volume Origin" within the 'System Preparation' tools pane.
4. Participant preparation in motion capture software program
   1. Within the main program window, open the 'Resources' pane. Click the "Subject" tab. Click the "Create a new Subject from a Labeling Skeleton" button. Select a labeling template from a list of template files provided.
   2. In the 'properties' window, enter the subject's name and values of body mass (kg), body height (mm), inter-ASIS distance (mm), left and right leg length (mm), left and right knee width (mm), and left and right ankle width (mm). In the 'Subject Resources' pane, right-click the subject name and click "Save Subject".
5. Record a static body calibration pose
   1. Ask the participant to stand motionless in the middle of the calibrated room with feet shoulder width apart while extending the upper extremities laterally so that all reflective markers on the body are well exposed to cameras.
   2. In the main program window, open the Tools pane. Click the "Subject Preparation" tab. In the Subject Capture section, click "Start" to record a 3-s motion trial to be the static calibration trial.
6. Procedure of determining functional joint centers
   1. Functional hip joint center
      1. Ask the participant to stand with one leg, and fully extend the other leg slightly forward. Instruct the participant to move the extended leg around the hip joint in the following sequence: move anteriorly and return to neutral, move anterior-laterally and return to neutral, move laterally and return to neutral, move posterior-laterally and return to neutral, move posteriorly and return to neutral, and a circumduction motion.
      2. Within the main program window, open the 'Tools' pane, click the "Capture" tab. In the Capture section, click "Start" to record a motion trial for each functional hip motion.
   2. Functional knee joint center
      1. Ask the participant to stand with one leg and maintain a 30° hip hyper-extension of the other leg. Instruct the participant to perform a 45° knee flexion with the non-weight bearing leg for 5 times.
      2. In the 'Capture' section of the 'Tools' pane, click "Start" to record a motion trial for each functional knee motion.
         NOTE: For details of the functional joint procedure, please see Schwarz, et al.²²
7. Motion capture of drop-landing movements
   1. Randomize the order of using three different drop-landing heights (26 cm, 39 cm, and 52 cm)¹⁴.
   2. Place the height adjusted wood box with a top surface area of 50 x 50 cm² on the floor covered by a rubber mat. The wood box is 11 cm from edges of the force plates. Ask the participant to stand on the box surface.
   3. Instruct the participant to extend their dominant foot directly in front of the box and shift their weight forward, and step off from the box. Ask the participant to land with both legs on the ground at the same time with each foot on a separate force plate.
   4. Ask the participant to remain standing until the motion capture of the trial is completed. Repeat the motion capture three times to collect three motion trials for each height.
8. Motion capture data processing
   1. Open a motion capture software program. Within the main window of the program, go to the 'Communications' pane. Click the "Data Management" tab. Select one of the recorded motion trials and open it in the program.
   2. In the 'Tools' pane, click the "Pipeline" tab. From the 'Current Pipeline' list, select the "Reconstruct" pipeline. Click the "Run" button to start the reconstruction process to obtain three-dimension (3D) trajectories of reflective markers.
   3. In the 'Tools' pane, click the "Label/Edit" tab. In the 'Manual Labeling' section, select individual marker names and label the corresponding 3D trajectories. Click the "Save" button of the toolbar when labeling is completed.
   4. In the 'Tools' pane, click the "Pipeline" tab. In the 'Available Operations' section, select "File Export". Double-click "Export C3D pipeline". Click the "Run" button to export the processed motion trial to a file in a coordinate three-dimension (C3D) format.
9. Biomechanical analysis of motion capture data
   1. Open a biomechanics software program to further process the motion capture data. From the top menu, click "File" and click the "Open/Add" button. Select the raw C3D files to import into the biomechanics software program.
   2. From the top menu, click "Model". Click "Create (Add Static Calibration File)". From the sub menu, select "Hybrid Model from C3DFile". Select and open the static calibration C3D file.
   3. From the top menu, click "Model". From the drop-down list, click "Apply Model Template". Select and open a model template file. Click the "Models" tab on the tool bar. Click the "Subject Data / Metrics" tab. Within the 'Subject Data' window, modify the values of 'Mass' and 'Height' to make the model subject-specific.
   4. Click the "Models" tab on the tool bar. Click the "Model Builder Advanced Post Processing" button of the top menu bar. In the pop-up window of the "Model Builder Advanced Post Processing", click the "Functional Joints" tab. Select "Add Motion File From Workspace".
   5. Select the functional joint center C3D files. Highlight an imported functional joint file. Highlight a functional joint matching the file. Use the "Set START Frame to Current Frame" and "Set END Frame to Current Frame" to select the appropriate portions of the motion trial. Click the "Compute Checked Landmarks" button. Repeat this process to calculate other functional joint centers to refine the skeletal model.
   6. Click the "Model" button on the top menu bar. Select "Assign Model to Motion Files". In the pop-up window of the "Assign Models to Motion Data", apply the subject-specific skeletal model to all motion trials.
   7. Click the "Pipeline" button of the tool bar. In the pop-up window of the "Pipeline Workshop", click the "Open Pipeline" button. Select the "Filtering Targets Pipeline". Click the "Execute Pipeline" button to perform a fourth-order low-pass Butterworth filter with cutoff frequency of 10 Hz on 3D trajectories of motion capture trials.
   8. Click the "Pipeline" button of the tool bar. In the pop-up window of the "Pipeline Workshop", click the "Open Pipeline" button. Select the "Filtering Forces Pipeline". Click the "Execute Pipeline" button to perform a fourth-order low-pass Butterworth filter with cutoff frequency of 60 Hz on ground reaction forces of motion capture trials.
   9. Click the "Settings" button of the top menu bar. Place check marks next to "Use Processed Analogs for Ground Reaction Force Calculations" and "Use Processed Targets for Model/Segment/LinkModelBased Items".
   10. Click the "Pipeline" button of the tool bar. In the pop-up window of the "Pipeline Workshop", click the "Open Pipeline" button. Select the "Model Based Calculation" pipeline. Click the "Execute Pipeline" button to perform calculations of lower-body joint kinematics and kinetics.
   11. Click the "Pipeline" button of the tool bar. In the pop-up window of the "Pipeline Workshop", click the "Open Pipeline" button. Select the "Export C3D Coordinates" pipeline. Click the "Execute Pipeline" button to export the processed 3D coordinates of lower-body visual markers in a C3D file.
   12. Click the "Pipeline" button of the tool bar. In the pop-up window of the "Pipeline Workshop", click the "Open Pipeline" button. Select the "Export Ground Reaction Forces" pipeline. Click the "Execute Pipeline" button to export the processed 3D ground reaction forces in a binary file (file extension: MAT).
       NOTE: To preserve the high impact peaks during landings, a cutoff frequency of 60 Hz is used to filter the raw ground reaction force data²³.
10. Preparing motion capture data for computer simulations
    1. Open a computer programming software. Import the filtered C3D data file and the MAT data file.
    2. Export a text file containing lower-body joint center coordinates. Convert the C3D data file and the MAT data file into text files (file extension: slf) for use by a multibody dynamic simulation program.

4. Subject Specific Modeling Procedure

1. Creating lower-body skeletal model
   1. Open the multibody dynamic simulation software program with the human body modeling plug-in installed. During this process, the human body modeling plug-in module is automatically opened. Within the splash screen, double-click the "New Model" icon to open the model building control panel.
   2. Within the main modeling panel, in the "Anthropometric Database Library" section, choose the generic body (GeBOD) from the drop-down list. Within the main modeling panel, specify body mass (kg), body height (mm), gender, and age (months).
   3. Within the main modeling panel, in the "Body Configuration" section, click the "Lower Body" radio button. From the "Units" drop-down list, select "Millimeter Kilogram Newton". Within the main modeling panel, click the "Apply" button in the "Create Body Measurement Table" section to accept the body measurements. Continue to click the "Apply" button in the "Create Human Segments" section to create a lower-body skeletal base model.
      NOTE: This model is scaled based on the individual's height, mass, age, and gender. The model consists of seven segments: a pelvis, two thighs, two shanks, and two feet (Figure 1). All of the segments are modeled as rigid bodies.
2. Modeling lower-body joints
   1. Within the main modeling panel, from the main menu drop-down list, select "Joints" to open the joint configuration panel.
   2. Within the joint configuration panel, in the "JOINT ROTATION ELEMENTS" section, click the button next to the "Prepare Model with Recording Joints". In the "Spring Dampers and Joint Limits Properties" section, enter the following parameters: Nominal Joint Stiffness of 1 Nmm/°, Nominal Joint Damping of 0.1 Nmm∙s/°, Joint Stop Stiffness of 3.38E7 Nmm/°. Continue to select "Left Leg" and "Right Leg" by checking the radio buttons next to the names. Click the "Apply" button to accept the joint configurations.
   3. Within the main modeling panel, from the drop-down list of the main menu, select "Workflow". From the drop-down list of the sub-menu, select "Gait" and "Calibrate". In the "Joint Center Data" section, enter the participant's lower-body joint center file.
   4. Click the "Load" button to import the data to modify the locations of joint centers. In the "Load Static Trial" section, enter the static calibration motion capture trial (in slf file format, generation described in steps 3.8-3.10). Click the "Load" button to import the file to parameterize the lower-body skeletal model.
      NOTE: By default, the hip joints are configured as spherical joints with three degrees of freedom, knee joints are configured as revolute joints with one degree of freedom, and ankle joints are configured as universal joints with two degrees of freedom.
3. Modeling skeletal muscles
   1. Within the main modeling panel, from the drop-down list of the main menu, select "Soft Tissues". From the drop-down list of the sub-menu, select "Create Base Tissue Set". In the "MUSCLE CONTRACTILE ELEMENTS" section, click "Prepare Model with Recording Muscle Elements".
   2. In the "GLOBAL RECORDING ELEMENT MUSCLE PROPERTIES" section, click the radio button of the "Updated 45 Muscle Set".
   3. In the "GLOBAL RECORDING ELEMENT MUSCLE PROPERTIES" section, accept the following default settings for muscle properties: Passive Stiffness of 0.4448 N/mm, Passive Damping of 1.75 E-2 Ns/mm, Muscle Resting Load of 0.4448 N. Check the radio buttons of "Left Leg" and "Right Leg" for muscle assignments. Click the "Apply" button to accept the configurations.
      NOTE: The 45 Leg Muscle set includes the following muscles: Adductor Brevis, Adductor Longus, Adductor Magnus (three groups), Biceps Femoris Long Head, Biceps Femoris Short Head, Extensor Digitorum, Extensor Hallucis, Flexor Digitorum, Flexor Hallucis, Gastrocnemius, Gemellus, Gluteus Maximus (three groups), Gluteus Medias (three groups), Gluteus Minimis (three groups), Gracilis, Hamstring, Iliacus, Lateral Gastrocnemius, Medial Gastrocnemius, Pectineus, Peroneus Brevis, Peroneus Longus, Peroneus Tertius, Piriformis, Psoas, Quadriceps Femoris, Rectus Femoris, Sartorius, Semimembranosus, Semitendinosus, Soleus, Tensor Fasciae Latae, Tibialis Anterior, Tibialis Posterior, Vastus Intermedius, Vastus Lateralis, Vastus Medialis.

5. Multi-body Dynamics Simulations

1. Performing inverse kinematic simulation
   1. Within the main modeling panel, from the drop-down list of the main menu, select "Workflow". From the drop-down list of the sub-menu, select "Gait" and "Trial". In the "Dynamic Trial Data" section, enter the file name of a dynamic motion capture trial (in slf file format) and click the "Load" button to import the data. Continue to enter the corresponding ground reaction force data file (in slf file format) and click the "Load" button to import the data.
   2. Within the main modeling panel, from the drop-down list of the main menu, select "_Analyze". Run the Reparameterize Analysis to adjust the model posture to match the posture in the beginning of the dynamic trial.
   3. Open the Simulation panel. Disable the effects of gravity and ground reaction forces. Choose the whole motion trial as the length of the simulation.
   4. Specify a simulation time step of 100 steps/s. Run an inverse kinematic simulation driven by the motion capture data. Save the inverse kinematic simulation analysis.
2. Creating a motion tracker agent
   1. Open the Motion Tracker Agent Creation panel. Accept the default tracker name: MA_Track.
   2. Set the Translational Stiffness and Rotational Stiffness as 10 N/mm and 1,000 Nmm/°, respectively. Set the Translational Damping and Rotational Damping as 10 Ns/mm and 1,000 Nmms/°, respectively. Set all translational and rotational degrees of freedom as Driven.
   3. Note. As only the lower-body model is used for the forward dynamic simulation, a motion tracker is necessary to account for the instability due to the lack of upper body motion.
3. Training leg muscles
   1. Open the Soft Tissues configuration panel. Choose Closed-loop Simple for the muscle model. Set the following parameters for the muscle model: Proportional Gain of 1.0E6, Integral Gain of 1.0E6, and Derivative Gain of 1.0E4.
   2. Select the inverse kinematic simulation analysis to be the target of the muscle training. Apply the muscle training.
4. Importing a flexible tibia
   1. Open the Flexible Body Import panel. Perform the Alignment Mapping with three known makers and their corresponding nodes on the surface of the flexible tibia.
   2. Choose the rigid tibia to be replaced by the flexible tibia. Select the MNF file representing the flexible tibia. Select the muscle attachment mapping file for reattaching leg muscles to the flexible tibia. Import the flexible tibia to the musculoskeletal model.
5. Performing forward dynamic simulation with the flexible tibia in place
   1. Open the Simulation panel. Enable the effects of gravity and ground reaction forces. Disable the effects of motion agents.
   2. Choose to run the simulation for the length of the whole motion trial. Set simulation time step of 100 steps/s. Run a forward dynamic simulation driven by trained muscles. Save the forward dynamic analysis.

6. Creating a Flexible Tibia Model

1. Creating a 3D surface mesh model
   1. Open an image processing program. Import CT slices in the DICOM format. Create a mask using the region growing method to separate bone tissue from the surrounding soft tissues.
   2. Search for CT slices where the tibia and fibula are connected. Separate the tibia and fibula by erasing the mask along the conjunction of the two bones.
   3. Create a second mask using the region growing method to only include the tibia bone. Go through the CT slices to uncover cavities existing in the tibia mask. Fill the cavities in the mask. Create a 3D tibia object based on the tibia mask. Export the 3D tibia object as a file in drawing interchange format (DXF).
2. Creating a finite element tibia model
   1. Open a FE analysis software program. Import the 3D tibia model file with the DXF extension.
   2. Perform the Sweep command to remove duplicated elements and nodes. Perform the Volume Mesh command to create a FE tibia model with hexagonal elements of 3 mm x 3 mm x 3 mm. Assign the following material properties to all elements: Young's modulus of 17 GPa, Poisson's ratio of 0.3, and density of 1.9E-6 Kg/cm³.
      NOTE: Material properties are assigned to each element with the assumption that bone tissue is isotropic within the ranges of strain experienced by bone during dynamic motions^24,25,26.
3. Creating a flexible tibia model
   1. Within the main control panel, click the "Geometry & Mesh" tab. Select "Geometry & Mesh". In the "Geometry & Mesh" pop-up window, in the "Mesh" section, click "Add nodes" to create two new nodes to represent the centers of the knee and ankle joints.
   2. In the main control panel, click the "Links" tab. Select RBE2's. In the RBE2's pop-up window, create link connections of type 2 rigid body element (RBE2) between the joint nodes and surface nodes on the knee and ankle surfaces.
   3. In the main control panel, click the "Boundary Conditions" tab. In the "Boundary Conditions" section, click the "New" button. Select "DOF_Set Nodes". In the "Boundary Condition Properties" pop-up window, create a boundary condition by assigning six degrees of freedom to each of the two RBE2 joint nodes.
   4. In the main control panel, click the "Loadcases" tab. In the "Loadcases" section, click "New", select "Adams Craig-Bampton"¹⁹. In the "Loadcase Properties" pop-up window, click "DOF-Set Nodes". Select the dofset_nodes created in the above step.
   5. In the main control panel, click the "Jobs" tab. In the "Jobs" section, click "New". Select "Structural". In the "Job Properties" pop-up window, select the loadcase created in the previous step. Click the "Job Results" button. In the "Results" pop-up window, select "Stress" and "Strain". Also select "Kilogram" for Mass, "Newton" for Force, "Millimeter" for Length, and "Second" for Time. Click the "Run" button.
   6. In the "Run Job" pop-up window, click the "Submit" button to submit the job for a FE simulation and to create the modal neutral file (MNF) of the tibia¹⁶.

7. Strain Data Analyses

1. Exporting bone strain data
   1. Open the post processor of the multibody simulation program. Load the Durability plug-in program.
   2. Open the simulation with the flexible tibia by clicking the simulation name. Export the maximum and minimum principal strains and maximum shear strain of the nodes representing the antero-medial aspect of the mid-tibial diaphysis.
2. Processing raw strain data
   1. Open a computer programming software for data processing. Import raw strain data. Apply a fourth order low-pass Butterworth filter to the raw data with a cutoff frequency of 15 Hz.

Representative Results

A healthy Caucasian male (19 years, height 1,800 mm, mass 80 kg) volunteered for the study. Prior to data collection, the subject reviewed and signed the consent form approved by the University Institutional Review Board before participating in the study. The experiment was conducted under the Helsinki Declaration. The experiment was performed based on the following protocol.

In order to verify the accuracy of the forward dynamic simulation, lower-body joint angles from the simulation were compared to the corresponding joint angles measured from the motion capture data processed by a biomechanics analysis program. A statistical analysis software was used to calculate cross-correlation coefficients of the comparisons. The cross-correlation calculation allowed 10 lags in both positive and negative directions. Each lag corresponded to a one time step in the forward dynamic simulation (0.01 s). The maximum cross-correlation coefficients were identified.

Visual inspection of Figure 2, Figure 3,and Figure 4 demonstrates the similarities between the joint angles produced with the experimental data and with the simulation data. Strong cross-correlation coefficients were found between the experimental and simulation joint angles at zero lag (Table 1).

The peak strains at the antero-medial region of the mid-tibial shaft during landing from three different heights are presented in Table 2. Among the three landing heights, the 52 cm landing condition demonstrated the largest peak maximum principal, peak minimum principal, and peak maximum shear strains. In addition, it was observed that, as the drop height increased, the peak maximum principal strains increased.

Figure 1: Subject-specific musculoskeletal model created in the present study. This lower body musculoskeletal model includes six rigid segments (pelvis, left and right femurs, left tibia, and left and right feet) and one flexible tibia (right tibia). 90 leg muscles are attached to the model. For visualization purpose, each muscle is represented by a coral color line. Joint centers are represented by light blue balls for right lower body and purple balls for left lower body. Please click here to view a larger version of this figure.

Figure 2: Joint angle comparisons (in degrees) between experimental motion capture data and simulation data for drop-landing from 26 cm height. Solid lines represent joint angles computed with experimental motion capture data. Dotted lines represent joint angles produced by multibody dynamic simulation data. Vertical lines represent moments of impact. Please click here to view a larger version of this figure.

Figure 3: Joint angle comparisons (in degrees) between experimental motion capture data and simulation data for drop-landing from 39 cm height. Solid lines represent joint angles computed with experimental motion capture data. Dotted lines represent joint angles produced by multibody dynamic simulation data. Vertical lines represent moments of impact. Please click here to view a larger version of this figure.

Figure 4: Joint angle comparisons (in degrees) between experimental motion capture data and simulation data for drop-landing from 52 cm height. Solid lines represent joint angles computed with experimental motion capture data. Dotted lines represent joint angles produced by multibody dynamic simulation data. Vertical lines represent moments of impact. Please click here to view a larger version of this figure.

	Droplanding Heights
	26 cm		39 cm		52 cm
Lower-body Joints	Cross-correlation Coefficient	Lag	Cross-correlation Coefficient	Lag	Cross-correlation Coefficient	Lag
Ankle	0.998	0	0.998	0	0.999	0
Knee	1	0	1	0	1	0
Hip	0.999	0	1	0	1	0

Table 1: Cross-correlation coefficients and lags from comparisons between joint angles produced based on motion capture data and joint angles produced from simulation data. One trial at each height was used for the comparisons. Zero lag indicates no difference in time when joint angles were produced between the two approaches.

	Droplanding Heights
Bone Strain (µstrain)	26 cm	39 cm	52 cm
Maximum Principal	1160	1270	1410
Minimum Principal	-659	-598	-867
Maximum Shear	893	870	1140

Table 2: Tibia bone strains at the antero-medial aspect of the mid-tibial shaft during drop-landing from three different heights. Maximum principal, minimum principal, and maximum shear strains are presented.

Discussion

The purpose of this study was to develop a non-invasive method to determine tibia deformation during high impact activities. Quantifying tibia strain due to impact loading will lead to a better understanding of tibia stress fracture. In this study, a subject-specific musculoskeletal model was developed, and computer simulations were run to duplicate the drop-landing movements performed in a laboratory setting. The effect of drop-landing height on tibial strain was examined. In this study, we observed that as the drop-landing height increased, so did the peak maximum principal strains. Also, among the three landing conditions, the 52-cm condition resulted in the highest peak maximum principal, minimum principal, and maximum shear strains.

Limited in vivo data are available in the literature with regard to the effect of drop-landing on tibia strain. Milgrom et al., reported the maximum principal strain ranging from 896-1,007 µstrain during landings from three different heights (26, 39, 52 cm)¹⁴. Ekenman et al. reported an average strain of 2,128 µstrain during landing from a 45 cm height¹³. The maximum principal strain from the computer simulations were between 1,160-1,410 µstrain during landing from three different heights (26, 39, 52 cm), which were higher than those reported by Milgrom et al. but were lower than that reported by Ekenman et al.^13,14

The following reasons may contribute to the difference in strain between the current and previous studies. First, demographic differences exist between the subjects in this and previous studies. We used a physically active male subject. Ekenman's study involved a female subject¹³. Milgrom's study included both males and females and reported the average strains¹⁴. Second, footwear may play a role in differences in bone strain.Lanyon et al. studied the effect of footwear on tibial strains, they found that walking and running barefoot resulted in greater strains compared to wearing shoes¹². The current study used a barefoot landing protocol, the strain values calculated were greater than those by Milgrom et al. study, which used a landing protocol with standard athletic shoes¹⁴. Third, alterations in landing strategy may also influence the tibial strain. In the present study, it was possible that the subject might choose a strategy such as increasing trunk flexion to help reduce the impact when the drop-landing height increased. This strategy could help protect the tibia from large strains. Milgrom et al. also suggested a possible protective strategy used by his subjects¹⁴. Fourth, there could be a slight difference in locations where tibial strain was monitored. Our study examined the bone strain at the antero-medial aspect of the mid-tibial shaft. In Milgrom et al., strains were recorded from the medial region of the mid-tibial shaft¹⁴. The sagittal plane bending moment on the tibia during landing may result in high maximum principal strain in places near the anterior regions of the tibial shaft. Nonetheless, our strain results appear to be comparable to results from previous studies and fall in the strain range (400-2,200 µstrain) reported by those in vivo studies^10,13,14.

The tibial strain values obtained from this non-invasive approach are influenced by the accuracy of the musculoskeletal model. Cross-correlations were performed to examine the experimental joint angle data and computer simulation data during drop-landings. Strong correlation coefficients were found between the experimentally measured data and computer simulation data. This indicates that the subject-specific model developed in this study can reasonably replicate the drop-landing movements. In addition, the tibial strains reported in this study were well below 3,000 µstrain, which confirms the assumption derived from other studies that the tibia bone deformation is linear during drop-landings^14,15. Thus, with the calculated strain data being in the linear range and excellent replications of landing movement patterns, we concluded that the strain data obtained from this non-invasive approach were reasonably accurate. Furthermore, the current study only recruited one subject to examine bone strain during drop-landings. Future studies could examine whether there is a dose response relationship between drop-landing heights and tibia bone strains by using a large sample size.

The significance of this study is to introduce an innovative non-invasive method of measuring bone deformation. This non-invasive approach addresses the limitations associated with the conventional in vivo strain gauge measurement, which could not be applied to a large sample of human subjects. In addition, the current proposed method addresses limitations associated with a previously reported non-invasive method^16,17, which was impacted by using limited kinematic data to drive the simulations and was only suitable for studying low ground impact movements such as walking. As tibia stress fractures remain high in the athletic and military populations, it is critical to study the effect of high impact physical activities (e.g., running, jumping, and cutting) on tibial bone responses. The current innovative non-invasive approach appears to be a feasible solution for conducting these studies. This will shed light on developing adequate physical training protocols for athletes and military recruits to reduce tibia stress injuries. Furthermore, this innovative non-invasive method presents an opportunity to evaluate bone strains in other bones inaccessible with implemented gauges such as the femur and navicular.

Important issues related to this non-invasive bone strain measurement must be addressed here. Firstly, a generic lower-body musculoskeletal model is created based on the individual's age, gender, body mass, and body height by using the GeBOD database²⁷. Experimentally measured spatial locations of lower-body joint centers are used to refine the musculoskeletal model. Compared to the generic model, this subject-specific modeling approach presents a better musculoskeletal model of the individual's physical structure. Future studies could consider developing a full body musculoskeletal model for upper body movement during multibody dynamic simulations.

Secondly, there are 45 muscles assigned to each leg in the model. Origins and insertions of the muscles are anatomically determined²⁷. A simple closed-loop algorithm is used to manage individual muscle's force production. Specifically, the change of muscle length history during dynamic motion such as landing is recorded via the inverse kinematic simulation. When the forward dynamic simulation is run, a PID controller was assigned to each muscle and used to regulate the necessary muscle force for duplicating the muscle length history recorded earlier. This simple closed-loop algorithm produces excellent results in replicating joint kinematics. However, this approach does not account for neural coordination among muscles with similar functions and could not account for co-contractions from antagonists. Future works may consider using a Hill-based muscle model, which consists of an active contractile element (CE) and a passive elastic element (PE). The Hill-based model integrates the muscle's force-velocity and force-length relationships to produce tension. The calculated muscle force can then be compared to EMG data for validation.

Thirdly, a subject-specific tibia model is created from CT images to represent the true geometry of the tibia bone under investigation. While CT imaging is the primary method to obtain the true geometry of the tibia bone, other imaging techniques such as magnetic resonance imaging (MRI) can also be used to produce the subject-specific tibia model. Also, the current modeling protocol assumes the material property of the tibia to be isotropic. A generic density value of 1.9E-6 kg/cm³ and a single Young's modulus of 17 GPa are assigned to all tibial FE elements. Future studies may consider obtaining density values from all regions in the tibia. This can be done by introducing a calibrated phantom during the CT scan. Bone density can then be calculated based on CT's Hounsfield units. Young's modulus of the bone tissue can be further calculated based on density data. Assigning subject-specific material properties to the tibial FE model will yield more realistic bone strain results through simulations.

Fourthly, a modal FE analysis is used to compute bone strains. During this modal analysis, frequency responses are computed to match mechanical loadings (linear and angular forces) imposed to the knee and ankle joints. A flexible tibia represented by an MNF file is generated from the modal FE analysis. This flexible tibia is introduced to the subject-specific musculoskeletal model to replace the corresponding rigid tibia. During the subsequent forward dynamic simulation, deformation of the flexible tibia at each time step is quantified. Compared to the traditional FE analysis, which computes the mechanical responses of an FE object consisting of thousands of degrees of freedom (thousands of elements and nodes) at each time step of motion, this modal analysis approach deals with far less numbers of degrees of freedom within the frequency domain (e.g., 12 loading conditions from the knee and ankle joints). With the modal analysis approach, computation time is significantly reduced from multiple hours/days to less than 1 h for a typical simulation. Besides the benefit of consuming less computer time, modal analysis approach is ideal for computing small deformation (< 10%) experienced by stiff structures such as bone tissue.

Finally, the advantages of the current non-invasive approach over a previously reported method^16,17 must be addressed here. A) Our musculoskeletal model is refined to possess more accurate lower-body joint centers, which are produced through the functional joint assessment²². However, the previous method defines joint centers for the model based on the Plug-in Gait procedure²¹ with the help of using a limited number of visual markers. B) This model incorporates 45 muscles to each leg compared to only 12 muscles used in the previous model. Increasing the number of leg muscles in the musculoskeletal model would improve the quality of the simulation. C) During the inverse kinematic simulation, the musculoskeletal model is driven by a set of 34 visual markers placed on the lower body, which allows better duplication of the actual movement. In contrast, the previous approach only uses 16 markers to drive the same simulation, and this may introduce numerical errors to the simulation. D) During the forward dynamic simulation, the real ground impact forces are applied to this musculoskeletal model to simulate the movement. However, the previous method is not able to incorporate ground impact forces in the simulation. Without using real ground impact forces during forward dynamic simulations, the previous method is limited to study low impact activities. The above steps we take to improve the fidelity of the subject-specific musculoskeletal model appear to be successful for examining tibial deformation during human movements. The addition of incorporating true ground impact forces in simulations proves to be necessary to study bone strain during high ground impact activities.

In conclusion, in vivo tibia bone deformation is normally measured by the conventional stain gauge method. This approach is associated with limitations such as an invasive nature, fewer volunteers, small bone surface areas being analyzed, etc. A novel approach employed multibody dynamic simulations with modal FE analysis was proposed in this study to quantify tibia deformation during drop-landings. It is evident that this approach can address the limitations inherited from the conventional strain gauge measurement. In addition, as this approach benefits from using real experimental kinematic and kinetic data, as well as a subject-specific musculoskeletal model and flexible tibia to perform dynamic simulation and modal FE analysis, it represents a huge improvement in the research protocol over a previously reported method. Thus, this non-invasive approach utilizing subject-specific data for multibody dynamic simulations combined with modal FE analysis could become a promising tool to study tibial deformation during dynamic motion. Future research could employ this method to study bone strains during high impact activities for a large cohort to study injury mechanisms of bone stress fractures.

Materials

Name	Company	Catalog Number	Comments
CT Scanner	GE Medical System	N/A	Light Speed VCT. For performing tibia CT scan.
Motion Capture System	Vicon Inc	N/A	Vicon FX40 high speed cameras. For performing 3D motion capture.
Force plates	AMTI Inc	N/A	Collecting 3D ground reaction forces
Vicon Nexus	Vicon Inc	N/A	Motion capture software program. For processing visual marker trajectory data.
Visual 3D	C-Motion Inc	N/A	Biomechanics analysis software. For computing 3D kinematics and kinetics of human movements.
MATLAB	Mathworks Inc	N/A	Computer programming software. For performing raw data filtering, data conversion, and data processing.
ADAMS 2012	MSC Software Inc	N/A	Multibody dynamic computer simulation program.
LifeMOD	Lifemodeler Inc	N/A	A software Plug-in in ADAMS. For building human body musculo-skeletal models.
MIMICS 13	Materialise Inc	N/A	Image processing program. A 3D modeling tool to process imaging data. For creating 3D tibia model from CT scans.
MARC 2012	MSC Software Inc	N/A	Finite element analysis software. For performing volumn meshing, generating tibia FE model, and running modal FE analysis.
SPSS 19	IBM Inc	N/A	Statistical analysis software.