| Detailed description |
|
|
Hi:
I've been using OpenSim 4.1 to put together a cervical spine model using the Python (3.7) API. In total, this is a 24 degree-of-freedom model: three rotational degrees of freedom at each joint level representing flexion-extension, lateral bending, and axial twist.
I'll start with a brief description of my suspected bug before I dive into my lengthy debugging procedure: I have reason to believe that, during inverse dynamics, C0-C1 lateral bending seems to be getting confused with C2-C3 axial twist---at least as far as my intervertebral disc models go.
I've been modelling the intervertebral discs as "ExpressionBasedBushingForce"s with exponential functions from the literature. I've verified that these functions are working as expected for each functional spinal unit (FSU) consisting of two vertebrae and the intervertebral disc between them. I have created OpenSim models for each individual FSU from C7-T1 to the skull-C1, and they seem to be working great!
The problem occurs when I try Inverse Dynamics on the whole cervical spine with some sample kinematics from the lab. In particular, I noticed that there was an abnormally large axial rotation moment (~4 Nm) at C2-C3 which did not exist at the other vertebral levels (~0.5 Nm). The problem persisted even when I zeroed out gravity. Coincidentally, when I disabled the intervertebral disc at that level, this large moment artifact disappeared. Thinking I found the answer, I tried running ForceReporter Analysis to calculate the moments in the Bushing elements: they were approximately 0.5 Nm, well below this aberrant 4 Nm moment. I even ran the C2-C3 kinematics through the FSU model on its own and repeated the inverse dynamics. Again, the 4 Nm moment did not manifest, and the inverse dynamics moments were in-line with what I would expect. Finally, I calculated the expected moment in the disc by hand given the exponential functions and the kinematics; my answer matched what the Bushing element was calculating in the ForceReporter. I also tried zeroing the motion at C2-C3 in the kinematics file and rerunning inverse dynamics on the whole cervical spine model; the same 4 Nm moment appeared at C2-C3, despite no motion occurring at that joint.
After a lengthy debugging procedure where I systematically zeroed each degree-of-freedom in the kinematics file, I found that I can reliably reproduce the erroneous 4 Nm moment simply by having lateral bending at the skull-C1 and no other motion in the file. Moreover, when I take the skull-C1 lateral bending angles and plug them into the exponential functions at C2-C3, I recover the 4 Nm moment that Inverse Dynamics was calculating. I suspected it might be an error in the ExpressionBasedBushingForce, so I replaced these elements with ExpressionBasedCoordinateForce to no avail--it returns the same 4 Nm axial rotation moment at C2-C3.
I've checked over the model's construction several times now and I cannot seem to find any reason for this confusion.
Attached I have a bare minimum representation of the model (cspine.osim) with a folder containing the obj-files of the bones (bones.zip). I have also included the kinematics (kinematics.mot) which results in the 4 Nm axial twist moment at C2-C3 (kinematics_inverse_dynamics.sto). Additionally I have the same kinematics file with no lateral bending motion between the skull and C1 (kinematics_no_lat.mot) which results in more reasonable moment estimates at C2-C3 (inverse_dynamics_no_lat.sto). Finally, I've included a kinematics file with only lateral bending motion at C0-C1 (kinematics_only_lat.mot) which, despite no motion or gravity, results in that tremendous 4 Nm moment at C2-C3 (inverse_dynamics_only_lat.sto). All of these are organized in a zipped folder called bug_report.zip.
I hope this bug report is thorough enough to help! Let me know if there's anything I can do on my end to assist!
Thanks,
Jeff Barrett
PhD Candidate, Biomechanics
University of Waterloo | Department of Kinesiology
|
|
