Open Knee(s): A Free and Open Source Library of Specimen-Specific Models and Related Digital Assets for Finite Element Analysis of the Knee Joint

There is a growing interest in the use of virtual representations of the knee for musculoskeletal research and clinical decision making, and to generate digital evidence for design and regulation of implants. Accessibility to previously developed models and related digital assets can dramatically reduce barriers to entry to conduct simulation-based studies of the knee joint and therefore help accelerate scientific discovery and clinical innovations. Development of models for finite element analysis is a demanding process that is both time consuming and resource intensive. It necessitates expertise to transform raw data to reliable virtual representations. Modeling and simulation workflow has many processes such as image segmentation, surface geometry generation, mesh generation and finally, creation of a finite element representation with relevant loading and boundary conditions. The outcome of the workflow is not only the end-point knee model but also many other digital by-products. When all of these data, derivate assets, and tools are freely and openly accessible, researchers can bypass some or all the steps required to build models and focus on using them to address their research goals. With provenance to specimen-specific anatomical and mechanical data and traceability of digital assets throughout the whole lifecycle of the model, reproducibility and credibility of the modeling practice can be established. The objective of this study is to disseminate Open Knee(s), a cohort of eight knee models (and relevant digital assets) for finite element analysis, that are based on comprehensive specimen-specific imaging data. In addition, the models and by-products of modeling workflows are described along with model development strategies and tools. Passive flexion served as a test simulation case, demonstrating an end-user application. Potential roadmaps for reuse of Open Knee(s) are also discussed. Supplementary Information The online version contains supplementary material available at 10.1007/s10439-022-03074-0.


A1. Modeling and simulation workflow Model development Image Segmentation
Target outcome Volumetric reconstruction of the tissues of interest, specifically the definition of the tissue boundaries, as binary image volumes (.nii; NIfTI) 8 and raw (without smoothing) triangulated surface representations (.stl) 29 in the same coordinate system as imaging data. Tissues include bones femur, tibia, fibula, patella; cartilage femoral, tibial (medial & lateral), patellar; menisci medial & lateral, ligaments anterior/posterior cruciate, medial/lateral collateral, patellar; tendons quadriceps. Segmentation algorithms The segmentation approach was primarily manual, requiring the modeler to use common labeling tools, such as, brush, pencil, etc. to paint or fill in the tissue region of interest in image volume. Depending on the tissue, multiple image volumes were used to confirm tissue boundaries or assist the tissue volume generation. Several 3D Slicer tools were also used during the segmentation process to assist gross segmentation before manual touch up. These are described in the following.

GrowCut Segmentation
From 3D Slicer documentation 1 : automata. The algorithm works by using a set of user input scribbles for foreground and background. For N-class segmentation, the algorithm requires a set of scribbles In this workflow, GrowCut was used primarily for gross segmentation of large tissue volumes.

Label Map Smoothing
From 3D Slicer documentation 1 : aliasing algorithm followed by a Gaussian smoothing algorithm. The output is a smoothed In the smoothing parameters outlined below, Sigma (Gaussian smoothing parameter) was chosen based on image resolution. For example, for cartilage MRI with resolution = 0.35 x 0.35 x 0.7 mm, sigma was set to 0.7. It should be noted that label map smoothing was used to assist segmentation as an intermediate step and should not be confused with preparation and smoothing of raw geometry for meshing (see section on Geometry Generation).

Joint Smoothing
From 3D Slicer documentation 1 : er smoothly, like jigsaw puzzle pieces.
It should be noted that joint smoothing was used to assist segmentation as an intermediate step and should not be confused with preparation and smoothing of raw geometry for meshing (see section on Geometry Generation).
Tissue-Specific Segmentation Procedures Below are the guidelines for segmentation procedures based on tissue types. It is important to use the specified input MRI, and locate the boundary as described for each tissue. However, the segmentation procedures described are only suggestions, and one may choose to highlight the anatomy using whichever tool they find most effective and efficient for their manual workflow. Once all tissues were segmented, the label maps were saved as NIfTI files (.nii) 8 and converted to raw triangulated surfaces (.stl 29 , in units mm) (without any further smoothing procedures). Slicer provides output options to accommodate these formats.

Geometry Generation
Target outcome Geometric reconstruction of the tissues of interest as smooth and watertight triangulated surface representations (.stl) 29 obtained from and in the same coordinate system as raw geometry; ready for volumetric meshing. Input Raw triangulated surface representations of tissues of interest (without filtering and smoothing) in .stl 29 format in image coordinate system.
Surface processing procedures A 5-step procedure (LVTIT, see below for descriptions) was used to process raw triangulated surface meshes. The process included staged smoothing approaches interleaved with surface reconstruction and resampling. The procedures aimed to generate uniform and watertight surface meshes that were smooth, volume-preserving and visually maintaining the geometrical shape of the tissues.
Smoothing Algorithms This section provides brief descriptions of the MeshLab 7 algorithms that were used during processing of triangulated surface meshes.

Laplacian Smoothing [L]
For each vertex in the mesh, a new position is chosen based on average position with nearest vertex (as described in built-in documentation in MeshLab 7 is the number of times the process is repeated.

Surface Reconstruction: VCG [V]
From built-in documentation in MeshLab 7 : lgorithm that have been used for a long time inside the ISTI-Visual Computer Lab. It is mostly a variant of the Curless at al. e.g. a volumetric approach with some original weighting schemes, a different expansion rule, and another approach to hole filling be informed in relation to original image resolution from which raw surfaces are obtained; e.g., for cartilage images with a resolution of 0.35 x 0.35 x 0.7 mm, use a world unit of 0.35, 0.5, or 0.7 mm.

Taubin Smoothing [T]
From built-in documentation in MeshLab 7 : iteration. Based on Gabriel Taubin, A signal processing approach to fair surface design, Sig Tissue specific processing parameters The table below specifies the parameters used in the surface mesh processing for each tissue for the surface meshes that were used in the generation of provided models. These values were considered a starting point and the process was -MeshLab defaults were used. The outcome of surface mesh processing needed to be checked (visually and when possible, quantitatively) to ensure that there was no significant loss of geometric features or volumes. If necessary, surface repairing was performed, i.e. using other MeshLab 7 tools, to ensure a manifold and watertight surface; sometimes faces needed to be reoriented such that surface normals all pointed outward. The final surface representation of the tissue was (.stl 29 , in units mm). It should be noted that these parameters were developed using one specimen and needed some specimen specific adjustments for the remaining specimens. Due to manual application of these parameters, some decisions made by the modelers may not be documented (several modelers with varying experience were involved in this process). File names typically consisted of some sequence of the smoothing step to provide an indication of which steps were used for a given surface, however exact details of parameters were not documented as the decisions were made based on visual inspections of how the surfaces looked after each step and whether any changes were required to create water tight meshes. Oks003 surface meshes were created using the same parameters as the additional surface meshes provided later in the document.

Additional surface geometries
To facilitate mesh convergence studies, additional multiple surfaces for each component were generated. Care was taken to use the following parameters consistently. An effort was made to document any deviations. However, it should be noted that the stochastic nature of software used may give different results when using the same parameters. For specimen oks003, additional surface meshes were created from the smoothed surface meshes created for the model.

Volume Mesh Generation
Target outcome Mesh definitions (nodes, elements, surfaces; node, element, face sets), template constitutive models (rigid bones; deformable other tissue). Outputs included tetrahedral volume meshes (cartilage, menisci, ligaments, tendon) in binary format (.med) 24 with node, face, element sets to facilitate model assembly, material property definitions, and assignment of tissue interactions.
Software requirements SALOME. SALOME 26 is an open-source software that provides a generic platform for pre-(cad, meshing) and post-processing for numerical simulation (free and open source LGPL license, see http://www.salome-platform.org/) SALOME 26 includes built-in scripting functionality using Python 32 , which was required to utilize Python scripts for mesh generation and annotation, and model assembly. SALOME 26 7.8.0 was used to support in-house Python 32 scripts.
StlToMed.py. In-house Python 32 script to generate meshes (surface or volume) and node, element, face sets using an XML based input file that points to .stl 29 files, e.g., tissue geometries, -in Python installation; developed using SALOME 26  interactions between tissues were constraints or contact, respectively, between opposing regions of tissues. Determination of surfaces, e.g. face sets, and/or node sets for these constraints was based on geometric principles relating to the pair of meshes based on: Proximity find all the nodes on part 1 that are within a prescribed distance (related to the element size) of any node on part 2. Normals select faces on part 1 that have normal vectors that point toward the barycenter of part 2. This can also be limited by the proximity condition if desired. Contains one mesh contains the other, select elements that have normal vectors that point inward or outward. All all the surfaces. This analysis was automated using the Python 32 scripts with SALOME 26 during mesh generation and annotation. The following table specifies all the tie constraints and contacts surfaces between tissues, and which geometrical principles should be used to relate them. Mesh generation and annotation -in Python 32 and the previously described connectivity file resulted in meshes of the tissue components along with node, element and surface definitions and node, element, face sets (for material property assignment, tie and contact definitions, and loading and boundary assignment). The script was automated, resulting in mesh files (in .med 24 format and in units of the input files) for each tissue component. The script starts with generating triangular surface meshes (for rigid objects bones) and tetrahedral (four node) volume meshes (for elastic objects other tissues) using geometric reconstructions of the tissue (.stl) 29 . Surface mesh discretization for rigid objects is identical to that of the input .stl 29 file. Surface discretization of elastic objects is identical to that of the input .stl 29 26 interface to NETGEN 3D mesh 19  Optimize the algorithm modifies initially created mesh in order to improve quality of elements. Optimization process is rather time consuming comparing to creation of initial mesh.

On
The script automatically generated an element set including all elements within the tissue, which was later used to assign material properties. Node and surface regions for tie and contact definitions were .0 was used to scale the approximate characteristic length of the surface mesh to determine the proximity threshold. For was not used.
Meshes of tissues were visually inspected in SALOME 26 . This inspection facilitated confirmation of appropriate definitions of node and surface regions to prescribe tissue connectivity, e.g., ligament insertion and origin sites, and contact, e.g. between cartilage. If necessary these sets were interactively edited in SALOME 26 by adding/removing nodes, faces, etc. from groups.

Template Model Generation
Target outcome Template model for finite element analysis in FEBio 11 format (.feb 13 , XML 20 based text file) including mesh definitions (nodes, elements, surfaces; node, element, face sets), template constitutive models (rigid bones; deformable other tissue), template interactions between tissue (tie constraints, contact), and template loading and boundary conditions (for rigid objects); created from geometric reconstruction of tissues.
Software requirements SALOME. SALOME 26 is an open-source software that provides a generic platform for pre-(cad, meshing) and post-processing for numerical simulation (free and open source LGPL license, see http://www.salome-platform.org/). SALOME 26 includes built-in scripting functionality using Python 32 , which was required to utilize Python scripts for mesh generation and annotation, and model assembly. SALOME 7.8.0 26 was used to support in-house Python 32 scripts.
MedToFebio.py. In-house script to generate a template FEBio 13 model using XML 20  Input tetrahedral volume meshes (cartilage, menisci, ligaments, tendon) in binary format (.med) 24 with node, face, element sets to facilitate model assembly, material property definitions, and assignment of tissue interactions.

Model Assembly
26 built-in Python 32 , the previously described connectivity file, and mesh files (outcome of mesh generation and annotation step) resulted in a template FEBio 13 model file (.feb 13 , version 2.5) for finite element analysis. The script is automated and resulted in the following contents written in the model file: Node definitions list of the node coordinates for each tissue Element definitions list of node connections to define the elements for each tissue Node sets, face sets, element sets, surface definitions groups of nodes, elements, faces labeled to be used in tie and contact pairs Material definitions o For rigid bodies (bones) only a place holder density is assigned (1e-9 tonnes/mm 3 consistent with spatial units of mm) o For elastic bodies (other tissues) a place holder density is assigned as 1e-9 tonnes/mm 3 and a Neo-Hookean material is defined using FEBio material type Mooney-Rivlin (uncoupled) 12 and setting c1 = 0.1 MPa; c2 = 0 MPa; k = 100 MPa (consistent with spatial units of mm) as a place holder.
Surface pairs to define connections (ties) or interactions between the surfaces, specifies he following criteria: o The surface on a rigid body should be the master surface.
o The larger of the two surfaces should act as the master surface.
o If the surfaces are of comparable size, the surface on the stiffer body should act as the master surface.
o If the surfaces are of comparable size and stiffness, the surface with the coarser mesh should act as the master surface.

Loading and boundary conditions
o Assignment of nodes on deformable bodies, e.g. ligament insertion and origin sites, to relevant rigid body as boundary conditions to tie the nodes to rigid bodies.
o Setting of six degrees-of-freedom kinematics (translations and rotations) of all rigid bodies to fixed as a place holder. Note that FEBio assigns these boundary conditions at the center of mass of the rigid body. 13 Load Data define the load curves that were Output requests

Module
Step and control parameters selected parameters to set simulation configuration for implicit static analysis using quasi-Newton incrementation with output at desired o time_steps = steps (default for steps is 10)

Model Customization
Target outcome Updated model for finite element analysis in FEBio 13 format (.feb 13 , XML 20 based text file) with modified tissue-specific constitutive models, additional joint stabilizers and convenience structures, implementation of in situ strain, implementation of anatomically based coordinate system and kinematic chains, modified loading and boundary conditions and, output requests to quantify joint kinematics-kinetics and tissue mechanics relevant to simulation of passive knee flexion ; generated from template model file (.feb) 13 . AnatomicalLandmarks_p3.py. In house script to utilize anatomical landmarks (described in an XML file containing the desired model properties) to calculate and add anatomical points and definitions to the model. To be used with Python 32 , source code available https://simtk.org/plugins/datashare/?group_id=485# (OKS_model) DOI pending.

Software requirements
Input Template (or customized) model in FEBio 11 format (.feb) 13 and tissue constitutive models and parameters relevant to knee mechanics (described below) as a text based input file.

Constitutive Models and Parameters
Material properties and the constitutive models were adapted from literature as described below.
Bone All bones (femur, tibia, fibula, patella) were assumed to be rigid bodies. Bones have a much higher stiffness than the other knee tissues. Rigid body assumption simplifies computation, therefore decreasing the computational cost and facilitates definition of joint kinematics and/or kinetics as loading and boundary conditions. A density of 1e-9 tonnes/mm 3 (identical to water; consistent with spatial units of mm). It should be noted that density assignment did not have importance on static simulations without the action of gravity. It is provided as a place holder.
Cartilage Cartilage was modeled as a nearly incompressible Neo-Hookean material defined by FEBio 13 setting C2 parameter of FEBio material type Mooney-Rivlin (uncoupled) 12 . This is a fairly would be adequate and computationally less challenging for joint level simulations while providing an opportunity to understand local mechanical environment on and within the cartilage. Cartilage constitutive model and coefficients were similar to previous modeling study 10 , which reported an -Hookean material coefficients are noted below.
1e-9 2.54 0 100 All units are consistent with spatial units of mm.
* Density assignment is a place holder; it did not have importance on static simulations without the action of gravity.

Ligaments and Tendons
Ligaments and tendons were modeled as nearly incompressible, transversely isotropic, hyperelastic material with a Mooney-Rivlin ground substance (Neo-Hookean by setting C2 = 0). This type of representation accommodates tensile dominant behavior of the ligaments dictated by their fiber alignment across their longitudinal axis. The parameters were identical to a previous modeling study 21 , which fitted data from literature. These values are noted below. Meniscus Menisci were modeled as nearly incompressible, transversely isotropic, hyperelastic material with a Mooney-Rivlin ground substance (Neo-Hookean by setting C2 = 0). This material model is seemingly more complicated than transversely orthothropic linear elastic models in literature. Yet, it allowed a convenient way to represent the behavior of meniscus which is largely dictated by its circumferential stiffness (based on fiber alignment) with the capacity to sustain compressive loading. For convenience, the parameters were identical to those used in a previous modeling study of meniscus 28 and they are noted below. Constitutive modeling of the menisci required specification of fiber direction. In FEBio, fiber direction can be specified for each element individually, in the ElementData section 13 . The meniscus fibers were oriented to the oval fitted to the bounding box.

Components for Stabilization
The following components were not tissues of interest, however they were incorporated in the model primarily to stabilize patella and represent the extensor mechanism. Mechanical springs and kinematic joints were used to represent these anatomical structures.
Quadriceps Tendon Attachment The quadriceps tendon was attached to the femur proximally.
Since the quadriceps muscles are not being modeled, this attachment were represented by a spring. The use of a spring instead allowed a weakly constrained movement of the proximal quadriceps tendon, while still preventing some separation of patella from the femur during flexion. The spring adds a small force pulling the patella to remain in contact with the femur bone. A spring constant of 1 N/mm was chosen, but can be calibrated if needed. The spring attached to the imaginary rigid body at the top of the quadriceps tendon on one end, and the femur bone on the other. The insertion origin on both rigid bodies was chosen as the origin of the imaginary rigid body.
Patellofemoral Joint Ligaments The medial and lateral patellofemoral ligaments (MPFL, LPFL) were modeled as discrete elements and using the force-displacement curve were defined as a tension-only linear springs.
. The insertion points of the patellofemoral ligaments were located as follows: Coordinates of anatomical locations, e.g, femoral epicondyles and patellar regions, were located in the MRI if their identification proved to be difficult to find on the meshes.

Application of In Situ Strain
Prestrain formulation of FEBio PreStrain Plugin is described in detail in literature 18 ; from this literature: -free to the prestressed reference configuration is represented by Fp, which will be referred to as the prestrain gradient. The total elastic deformation gradient Fe is determined by the composited deformation To scale the in situ stretch as expected when reducing the step size , a load curve was defined for each prestrain ligament to set in situ stretch to 1 (no strain) at time = 0. For example, if the desired initial stretch value is 1.034, the load curve should be: <loadcurve id="1" type="linear"> <point>0,1.0</point> <point>1,1.034</point> </loadcurve> To use this load curve appropriately, the initial stretch defined in the material was changed to 1.0, e.g. <stretch lc="1">1.0</stretch> One average initial stretch was defined for all fibers in each ligament as:

Ligament
Initial Stretch  9 , who referred to previous knee models 21 and experimental data 14 .
# Initial strain set to zero due to lack of data, similar to Dhaher et al. 9 In the prestrain constraint, min_iters and max_iters were set to 0, essentially eliminating augmentation. Simulations solely rely on penalty based constraints to balance the higher possibility of convergence with potentially decreased constraint enforcement.
Since the prestrain update rule is chosen to eliminate distortion, the above values would only be used as an initial guess for fiber stretch. Due to changes in cross sectional area of the ligament, in order to maintain equilibrium, the solver updates the fiber stretches as needed. Due to this in previous studies 9,21 ), as the values are likely to change. Thus, one average value was given as the initial guess for fiber stretch for all regions in the ligament. After the solver determined the updated initial stretch values, they can be compared to the above literature values.

Locating Joint Axes
Tibiofemoral Joint The Grood and Suntay Joint Coordinate System 15 (JCS) was used to define this joint. This method requires a definition of a tibial coordinate system (xT,yT,zT), a femoral coordinate system (xF,yF,zF), and a floating axis (FTF). The following anatomical landmarks were located on the meshes of the femur and tibia: Flexion Axis (xF-axis): The cross product between the yF-axis and the zF-axis.
Tibiofemoral Floating Axis: The tibiofemoral floating axis (FTF-Axis) is defined as the cross product between the zT-axis and the xF-axis at any given joint position.
Patellofemoral Joint This joint motion is described in a method similar to JCS 15 , according to Bull et al. 5 This method requires a definition of a patella coordinate system (xP,yP,zP) , femoral coordinate system (same as that described above for the tibiofemoral joint), and a Floating axis (FPF ). The following anatomical landmarks were located on the mesh of the patella: Anterior-Posterior Axis (yP-axis): The cross product between xP-axis and zP-axis.
Patellofemoral Floating axis: The patellofemoral floating axis ( FPF-axis) is defined as the cross product between the zP-axis and the xF-axis.
Note:. We defined the superior-inferior axis as perpendicular to zp, and in line with the coronal image plane (y=0). This would not work if the knee is rotated in the MRI coordinate system. If this n to help -femoral joint constraints are all set free, this should not affect the model, as the joint has 6 degrees of freedom.
Creating Kinematic Joints in FEBio Recent versions of FEBio 13 do not allow independent prescription of rotational degrees of freedom for a rigid body. For this reason and to facilitate loading and boundary conditions relevant to anatomy of the knee, kinematic joint chains were defined for tibiofemoral and patellofemoral joints.
Tibiofemoral Joint: Three rigid cylindrical joints 13 were added to the constraints section to define tibiofemoral were defined: TFTO: located at the origin of the tibial (xT,yT,zT) coordinate system.
TFFO: located at the origin of the femoral (xF,yF,zF) coordinate system.
The cylindrical joints were defined as follows:

FTF-axis TFFO TFTO
All other parameters in the rigid cylindrical joints were set to FEBio defaults 13 . Each cylindrical joint also had a translational component to describe tibiofemoral joint translations 13 .
Patellofemoral Joint: Three rigid cylindrical joints 13 were added to the constraints section to define patellofemoral were defined: PFPO: located at the origin of the patella (xP,yP,zP) coordinate system.
PFFO: located at the origin of the femoral (xF,yF,zF) coordinate system.
The cylindrical joints were defined as follows: Note: Centers of mass of tibia, patella, femur rigid bodies were moved from world origin to their respective coordinate system origins, in order to allow for easier manipulation while debugging the model, and for easier calculations in post-processing.

Loading and boundary conditions
The simulation strategy was to use a one step solution. Load curves defined the in situ strain (prestrain) application from time 0 to 1, then hold the prestrain at its final value between time 1 and 2. The loading (flexion) occurs from time 1 to 2, and the load curve for the flexion is set to be zero from time 0 to 1, i.e., fixing the flexion angle during the prestrain step.
In Situ Strain Application In situ strains were applied based on previously described specifications (see above). Any loading and boundary conditions that are not specified below were set to FEBio defaults 13 .
Prestrain: as described in in situ strain application section (see above) Femur: All degrees of freedom (3 translations, 3 rotations) free.
Passive Flexion At the start of the passive flexion application, in situ strains should have already been applied as prescribed. Our interpretation of passive knee flexion is that the motion of the knee is guided by joint contact and connective tissue recruitment. Therefore no external loading was applied; the movement of the knee was unconstrained other than prescription of the flexion angle up to 90º. Any loading and boundary conditions that are not specified below were set to FEBio defaults 13 .
Quadriceps tendon slider joint: All degrees of freedom (1 translation) free.

Simulation output requests
By default, FEBio 13 provides two output files: .xplt, plotfile, a binary file that includes all default and requested field variables compatible with PostView 22 ; .log, logfile, a text file that reports the convergence history of the simulation along with convergence metrics including any other variables requested for output. In accordance with the overall customization of the model and to provide output metrics relevant to simulation case, output of the following metrics were requested: Fiber stretch fiber stretch during simulation in tissues with compatible constitutive models (ligaments, tendon, menisci).
Generalizable parameters contact gap, pressure, traction; displacement; reaction forces; stress; to be stored in plotfile.
Prestrain stretch actually the total fiber stretch, including the effects of the prestrain and deformation; to be stored in plotfile.
Fiber stretch fiber stretch due to deformation during simulation in tissues with compatible constitutive models; this differentiation is important to understand ligament deformations when in situ strain is implemented.
Rigid body data Request outputs kinematics and kinetics of all rigid bodies including imaginary rigid bodies used for joint coordinate systems; to be stored in plotfile and logfile. Information in logfile can be used to reconstruct joint kinematics.

Rigid joints data
Request outputs only kinetics but not kinematics as it is not implemented in FEBio 11 ; to be stored in plotfile and logfile. 13 for more details.

Simulation
Target outcome Solution of fully customized model through finite element analysis using FEBio 11 ; generating simulation results as binary and text output files (.xplt and .log, respectively).
Software requirements: FEBio. FEBio 11 is a nonlinear implicit finite element analysis framework designed specifically for analysis in biomechanics and biophysics (custom open source license; free for academic research use, licensing for commercial use is available, see http://www.febio.org). The latest version of FEBio 11 (version 2.9 at the time of preparation of this document) was used.
FEBio PreStrain Plugin. PreStrain 18 Plugin provides a general framework for representing prestrain in a finite element model using a prestrain gradient method. A version compatible with the latest version FEBio was used (at the time of preparation of this document, version 1.0 supporting FEBio 2.9 has been available).

Input Fully customized model in FEBio format (.feb 13 ).
Simulation Process Invoked FEBio 13 with the model file as input. For simulations that did not convergence, relaxation of convergence tolerances and utilization of alternative solution algorithms, contact formulations, etc. were employed (see Appendix A6).

Post-Processing
Target outcome Extraction and summary of knee kinematics and kinetics during passive flexion; processed using raw simulation results of fully customized model with FEBio 13 (.log file), supported by graphs as binary image files.
Input Solution of fully customized model through finite element analysis using FEBio 13 ; simulation results as binary and text output files (.xplt and .log, respectively).

Kinematics of tibiofemoral cylindrical joints: 3 rotations, 3 translations total
Kinematics of patellofemoral cylindrical joints: 3 rotations, 3 translations total Translation of tibia origin relative to femur origin in femoral coordinate system Translation of patella origin relative to femur origin in femoral coordinate system Constraint moment to maintain the knee joint at prescribed flexion Visualization PostView 22 was used to take snapshots of the model at different flexion angles, as obtained through simulation of passive flexion. PostView 22 can also be used to inspect tissue stress-strain distributions, export data, images, and animations.

A2. Experimental data processing
Registration for Specimen-Specific Calibration Target outcome Coordinate system transformation matrices between joint testing and imaging coordinate systems of bones and experimental anatomical landmarks transformed to model coordinate system in XML 20  Input Experimental probed points on registration markers and anatomical landmarks and coordinate system transformations (State.cfg); raw registration marker geometries from imaging; FEBio 13 model file.
Registration In the State.cfg file, one can find the coordinates of the probed points on registration markers (three on femur, three on tibia) and bone landmarks collected during mechanical testing 6 . A sphere was fit to probed points on each registration marker to obtain its center in the local bone motion tracking system coordinate system. Similarly, a sphere was fit to raw registration marker geometries obtained by segmentation of imaging data to obtain their centers in image coordinate system. For each cluster of registration markers on the bone, the transformation matrix was calculated between the local bone motion tracking system and image coordinate system using singular value decomposition between sphere centers. In following, anatomical landmarks on each bone, which are probed during mechanical testing, were transformed to image coordinate system to serve as the foundation to redefine model joint coordinate systems and cylindrical joint axes. The coordinate systems of tibia and femur were updated based on descriptions provided in the experiment documentation 16 (also see below). Patella coordinate system remained the same due to incomplete data on patella registration marker assembly. Customized Full Knee Model with Experiment Coordinate Systems Experimental landmarks were used to define anatomical joint coordinate system in the model. The tibiofemoral floating axis (FTF-axis) was defined as the cross product between the Tz-axis and the Fx-axis at any given joint position. The patellofemoral floating axis ( FPF-axis) was defined as the cross product between the Pz-axis and the Fx-axis. The in house script FebCustomization_p3.py needs to be run for this purpose. The model were therefore aligned with the experiment such that the axes as defined in the experiment are the same as the ones defined in the model (For script usage, see Appendix A8).
Input Experimental joint kinematics-kinetics data (.tdms), joint coordinate system offsets (State.cfg); coordinate system transformation matrices between joint testing and imaging, and experimental anatomical landmarks transformed to model coordinate system (.xml).
Processing of Passive Flexion Data Experimental passive flexion data file (.tdms) were processed. Kinematics data were 3 , which provide joint kinematics in an anatomical joint coordinate system (defined in the experiment based on cylindrical joints) 16 relative to a reference state (joint offsets given in State.cfg 4 ). Kinetics data were 3 , which provide joint kinetics in an anatomical tibia coordinate system as applied as loads on tibia 16 . The Python script extracted, processed, and stored the data: 1. Extracted data such that the data from 0° to maximum flexion was re-sampled at 5° increments, averaging data on each channel where flexion angle was within +/-0.5 °.
2. Added kinematics offsets (from State.cfg) to kinematics channels to report bone pose and orientation in an absolute fashion.
3. Transformed kinematics data to the convention used in the model, i.e., cylindrical joint translations and rotations, accommodating offsets at model reference state when reconstructing experiment coordinate systems in the model. 4. Transformed kinetics data to the convention used in the model, i.e., joint loading applied to femur in model coordinate system, which is registered and aligned to experiment coordinate system.
Thresholds for cropping and resampling of data may change depending on data quality, e.g., noise and errors that may become apparent during analysis. The content of experimental kinematics-kinetics data files are in right knee abstraction, which were managed during any coordinate system transformation.
Processing of Laxity Data Kinematics-kinetics data for laxity were split into files based on flexion angle (0º, 30º, 60º, 90º), dominant loading axis (anterior-posterior translation, internal-external rotation, varus-valgus), and direction of loading axis (positive, negative). Kinematics data were 3 , which provide joint kinematics in an anatomical joint coordinate system (defined in the experiment based on cylindrical joints) 16 relative to a reference state (joint offsets given in State.cfg 4 ). Kinetics data were extracted from 3 , which provide joint kinetics in an anatomical tibia coordinate system as applied as loads on tibia 16 . A Python script was developed for extraction, processing, and storage of laxity data: Extracted data by finding the indices of the data points (from Kinetics.JCS.Desired) where the force was held constant for all loading cases at all degrees of flexion where laxity data was collected (0 º, 30 º, 60 º, 90 º), taking the index of the final data point for each flat section of the data and ext Added kinematics offsets (from State.cfg) to kinematics channels to report bone pose and orientation in an absolute fashion.
Transformed kinematics data to the convention used in the model, i.e., cylindrical joint translations and rotations, accommodating offsets at model reference state when reconstructing experiment coordinate systems in the model.
Transformed kinetics data to the convention used in the model, i.e., joint loading applied to femur in model coordinate system, which is registered and aligned to experiment coordinate system.
Thresholds for cropping and resampling of data may change depending on data quality, e.g., noise and errors that may become apparent during analysis. The content of experimental kinematics-kinetics data files are in right knee abstraction, which were managed during any coordinate system transformation.
A3. Model customization for application of kinematics-kinetics (with or without experimental data) Target Outcome Customized full knee models in FEBio 13 format (.feb, XML 20 based text file) prepared for all simulation cases, and all experimental loading conditions. Models include converged meshes, confirmed material properties and calibrated in situ ligament strains, and for reproduction of experiments, loading and boundary conditions of joint testing registered and transformed to model coordinate system.
Customization for Experiment or simulated Loading Cases Python script updates knee model in FEBio 13 to replicate experimental or simulated conditions. Experimental kinetics were applied as external femur loads, and experiment flexion angle was prescribed to the extension-flexion joint. Tibia was fixed; femur and patella were free to move and all loads and boundary conditions were applied in one step. From time 0 to 1, in situ strain were applied while keeping flexion at 0°. From time 1 to 2, the loads and boundary conditions at the start of experiment were prescribed, i.e., the flexion angle was set and the loads in the remaining degrees of freedom were applied on femur. From time 2 to 3, the loads and boundary conditions of the experimental trial were applied until the end point of the experiment. Load curves for each degree of freedom (particularly the dominant loading) were defined based on experiment data points and simulation output were requested at each experiment point. The kinematics-kinetics trajectories of experiment were split to facilitate prescription of loading scenarios in simulations.

A4. Model calibration procedure Mesh Convergence
Target Outcome Geometric reconstruction of tissues of interest as smooth and watertight triangulated surface representations (.stl) 29 and finite element meshes (.med) 24 with several mesh densities including converged mesh densities. Compartmental models of tissues of interest in FEBio 13 format (.feb, XML 20 based text file) for mesh convergence simulations. Simulation results as binary and text output files (.xplt and .log, respectively) 13  Mesh Generation Procedures For each tissue geometry that was generated above, a mesh was created in Salome 26 . A new MED 24 file was generated where all node/element/face groups, which were already created for original tissue mesh, were transferred to the new MED 24 file. This was done in a scripted fashion, using in house Python 32 script transfer_med_groups.py by treating the groups as a binary field and applying that binary field to the new mesh to find the correspondence of nodes and elements in the new mesh. This correspondence was used to select nodes and elements to generate groups with the new node and element sets. A quality assurance check was completed here to ensure that all node/element/face groups are as expected. For example, a check was done for all ligament insertion origins that the insertion area covers the entire width of the ligament where it connects to the bone.
Model Generation Procedures Each of the tissues of interest was tested by creating a model with the boundary conditions outlined in the table below. The models were created for each mesh density, such that 4 models exist for each tissue, with the only difference being the density of the mesh of the tissue of interest. Model generation was performed according to the model development A1 (using in home scripts MedToFebio.py, FebCustomization_p3.py), and boundary conditions were updated manually in the FEBio model input files for each model.   Mesh Convergence Procedures Each of the models created above was run in FEBio, beginning with the coarsest mesh density. The simulation was repeated, each time with a finer mesh. The primary measured outputs were compared with those of the previous simulation, and when there was less than 5% difference (7.6% for medial tibial cartilage), convergence was assumed. The mesh density at which the output converged was then used in the final model. If no convergence is found within the specified mesh densities, finer/coarser mesh densities were created to continue testing until convergence was found.
Template Full Knee Model with Converged Meshes The in house scripts MedToFebio.py, and FebCustomization_p3.py were run including all of the selected knee tissues with converged mesh densities.

Confirmation of Material Properties
Target Outcome Compartmental models of tissues of interest in FEBio 13 format (.feb, XML 20 based text file) for simulations to confirm and modify material properties using converged meshes. Simulation results as binary and text output files (.xplt and .log, respectively) 13 and as summary of calibration process including target metric, model predictions as a function of material property and fit error (XML 20 based text file). Full knee model with confirmed and modified material properties in FEBio 13 format (.feb, XML 20 based text file). Tissues for which material properties were confirmed include cartilage femoral, tibial (medial & lateral), patellar; menisci medial & lateral, ligaments anterior/posterior cruciate, medial/lateral collateral, patellar; tendons quadriceps.
Input FEBio model file of the full knee with converged meshes, including all the tie/contact surfaces.
Compartmental Modeling of Structural Tissue Behavior In order to confirm and modify material properties, structural tissue response or gross material behavior were compared with literature. If, the tissue behavior fell outside the range of expected behavior, the material properties were adjusted until the mechanical response of the tissue is considered within the reported normal range. This process essentially confirmed each tissue material property used for model development with a second information resource. Details of the simulations (and/or analyses) performed are provided below. An in house script ModelReduction_rigid3.py, was used to reduce the full knee model to include only the desired components for each of the test simulation cases. An in house script StiffnessFromLog.py was used to determine linear stiffness from model results. For the tension test simulations in the fiber direction, a load curve was defined from 0 to X, where X is the desired total displacement. Then, in the boundary section the prescribed displacement was specified in the fiber direction. For example, given a fiber direction of [-0.101,-0.511,0.854] a unit displacement (X=1) : <rigid_body mat="4"> <prescribed bc="x" lc="3">-0.101</prescribed> <prescribed bc="y" lc="3">-0.511</prescribed> <prescribed bc="z" lc="3">0.854</prescribed> <fixed bc="Rx"/> <fixed bc="Ry"/> <fixed bc="Rz"/> </rigid_body> Python script StiffnessFromLog.py was used to extract rigid body reaction force displacement data from simulation output (.log). In following, force-displacement response of the tissue along its fiber direction was calculated. Linear stiffness of the tissue (k) was calculated by fitting a line (E) was calculated as (k x Lo / Ao), where Lo and Ao are the reference ligament length and cross sectional area, respectively. The reference length was calculated as the distance between the centers of the insertion and origin. The cross sectional area was calculated as the average cross sectional area between the insertion and origin. If the simulated tissue properties (k and/or E) fall outside two standard deviations of the expected behavior, the fiber modulus (C5) was scaled to bring the tissue response within the expected range.

Cartilage
Since we assume cartilage properties to be consistent for all cartilage tissues, we performed one indentation test to confirm if cartilage structural response is within what is reported in literature. For this purpose, we used indentation stiffness of lateral tibial cartilage, which was reported as 20.38±5.32 N/mm 30 . A model was generated to reproduce the experiment conditions 30 . The model included the tibia bone (TBB), lateral tibial cartilage (TBC-L), and a 1mm diameter indenter as a rigid body, which was in frictionless contact with the cartilage. The indenter was placed above the cartilage near the meniscus, as described in the study, and a load of 0.5 N was applied to force the indenter against the tibial cartilage. StiffnessFromLog.py was used to extract indenter force and displacement data from simulation results (.log) and the linear stiffness was calculated using the linear region (upper third) of the force-displacement curve. Material properties of cartilage (C1) were scaled as needed to get the range within two standard deviations of the reported stiffness value. Bulk modulus parameter (K) scaled accordingly.

Meniscus
Depending on the location of the sample (anterior, central, posterior) and the thickness of it, circumferential tensile modulus of medial meniscus was reported as 43.4±26.8 MPa to 141.2±56.7 MPa 17 . The fiber modulus (C5) of meniscus was scaled to match within two standard deviations of the reported modulus.

Customized Full Knee Model with Confirmed Material Properties
The in house script FebCustomization_p3.py needs to be run including all of the selected knee tissues with updated material properties.

Calibration of In Situ Ligament Strains
Target Outcome Full knee models with converged meshes, confirmed material properties, joint coordinate system defined to align with the experimental coordinate system, and loading and boundary conditions of experiments selected for calibration in FEBio 13 format (.feb, XML 20 based text file). Simulation results as binary and text output files (.xplt and .log, respectively) and as summary of calibration process including target metric, model predictions as a function of in situ ligament strains and fit error (XML based text file). Full knee model with calibrated in situ ligament strains in FEBio format (.feb, XML based text file). Tissues for which in situ ligament strains were calibrated include ligaments anterior/posterior cruciate, medial/lateral collateral. Input Template FEBio model file of the full knee (.feb 1 ) and model properties (.xml 20 ) files with converged meshes, confirmed material properties, and experiment coordinate systems; processed specimen-specific kinematics-kinetics data (.csv 27 ). Models experiment_to_model.py, was used to generate models representative of the loading and boundary conditions of selected laxity tests to calibrate in situ strains. Only joint laxity data at 0° flexion was used to modify in situ strains for anterior cruciate ligament (ACL), posterior cruciate ligament (PCL), medial collateral ligament (MCL), and lateral collateral ligament (LCL). The decision to use only these data was made to save on computational cost and time. All other loading scenarios include flexion of the joint prior to performing laxity testing, which can be costly, and often, convergence issues may arise. This way, the calibration could be performed quickly, and the models were unlikely to have any convergence issues.
Template model (.feb) and model properties (.xml) reflective of converged meshes, confirmed material properties, and experiment coordinate systems were the basis for customization of models for calibration. Modifications to the customization script and/or additional scripts were necessary to implement the loading scenarios Overall, application of loading and boundary conditions and output requests was similar to those described in A1 with exceptions noted in here. Tibia was fixed; femur and patella were free to move and all loads and boundary conditions were applied in one step. From time 0 to 1, in situ strain was applied while keeping flexion at 0°. From time 1 to 2, the loads and boundary conditions at the start of experiment were prescribed, i.e., the flexion angle was set and the loads in the remaining degrees of freedom were applied on femur to reflect the loading state of the joint at the start of testing. This step accounted for any offsets in bone configuration between imaging and the experiment and it should be done after the prestrain step as we do not want the in situ strain calibration to be dependent on the orientation of the knee in different experiment trials. From time 2 to 3, the loads and boundary conditions of the experiment were prescribed, i.e., the flexion angle was constant and the loads in the remaining degrees of freedom were applied on femur. Load curves for each degree of freedom (particularly the dominant loading) were defined based on experiment data points and simulation output were requested at each experiment point. A total of 4 models were generated:

A6. Model specific measures for convergence
Once the models were customized, some model specific adjustments had to be made for full convergence. Models were run in FeBio 2.9. If run in other versions, due to inherent software behavior, other measures might be needed to achieve convergence.
oks001 Stiffness for springs between medial meniscus and MCL was reduced from 1000 N/mm to 300 N/mm for full convergence.
oks002, oks004, oks008 Stiffness for springs between medial meniscus and MCL was reduced from 1000 N/mm to 3 N/mm for full convergence.
oks006 Adjusted ties for MCL-femur. Stiffness for springs between medial meniscus and MCL was reduced from 1000 N/mm to 3 N/mm, MCL prestrain reduced to 0.025 and, LCL prestrain reduced to 0.02.
oks007 Stiffness for springs between medial meniscus and MCL was reduced from 1000 N/mm to 3 N/mm and MCL prestrain updated to 0.025 for full convergence.
oks009 Stiffness for springs between medial meniscus and MCL was reduced from 1000 N/mm to 3 N/mm, updated MCL-Femur ties for full convergence.
Note check FeBio_custom.feb model file for inactive contact definitions for each model.

A7. Database folder structure
Various model related data are organized in subfolders under folder with the specimen name.

tdms_processing_oks.py
Python script used to process the tdms files from open knees data. Saves the data as a png graphs and csv files. plot_groups function can be used to save the csv, png files in intermediate steps (raw, partially processed, etc).

experiment_to_model.py
This script was designed to take processed experimental data (kinematics in JCS with model offsets already removed, kinetics as forces/moments applied to the femur defined in the tibia coordinate system), and the model that results from running the customization script, and, create a model that applies to experimental conditions to the model by first converting the kinetics data into the image coordinate system, and then applying those kinetics as external forces to the femur. Kinematics load curves were also added, so that the kinematics can be easily applied to the cylindrical joint by setting the translation/rotation of the joint to follow the load curve. The make_model function will do this for one model. NOTE -it will change the existing febio_file if name is not given make_model(modelprops_file, febio_file, kinetics_csv, kinematics_csv, name=none) To run it for several files at the same time, read_from_xml(exp_to_mod_xml) function should be used.

InSituOptimization.py
This script was designed to perform calibration of the ACL, PCL , LCL, MCL in situ strains. The input for this function is an xml file that points to the models used in the calibration, and the , https://simtk.org/plugins/datashare/?group_id=485# (OKS_model) DOI pending. Under the Files section: provide the path to the model properties file. The model properties file should include the added anatomical landmarks from running the customization script, as the joint axes which are added are used to process the model kinematics. Under the Options section: rms_error type can be defined as "loading" or "all". loading = minimize rms error in the loading direction only (ex for ACL find in situ strain to only minimize mismatch on anterior translation axis) all= minimize rms error on all 6 dof (ax for ACL in situ strain to minimize mismatch on all translation and rotation axes). opt type can be defined as "single" or "multi" . single = single variable optimization, optimize in situ strain for one ligament at a time. multi = multi variable optimization, optimize in situs train for all ligament simultaneously (this will take much longer) Under the Ligaments section: for each of the ligaments, provide the febio file (created using experiment_to_model.py). For example, for the ACL, this would be a model generated to mimic the 0 degree flexion anterior loading test. (doesn't need to be 0 degrees, we just decided 0 because the models run the fastest) provide the experiment kinematics file -needs to match with whatever file was used for the febio file. For example, if ACL 0 degree anterior laxity, then this should be the kinematics file for that trial in situ strain is the initial guess for in situ strain. This shouldn't matter in theory -eventually it should find the solution. However, if the optimization has already run and a small change is made and need to run it again to be sure it doesn't affect results, etc. a known solution can be put as the initial guess to save a lot of computation time. The order in which these ligaments appear in the xml file is the order in which they will be optimized. Running this script will create summary files named Optimization<num>.txt, where num is the counter for which round of optimization it is on. The script will continue to loop on all the ligaments until convergence is achieved. Essentially create febio models for each laxity test (using experiment_to_model.py) and put them each in their own folder. Then create the InSituOpt.xml file pointing to each of these files, and other relevant files. Then run from the terminal as: python InSituOptimization.py /path/to/InSituOpt.xml febioCommand febioCommand is the path to the febio execulatble, on linux it can be found in '/FEBio2.8.2/bin/febio2.lnx64' If always running from the same computer, the febio command can be hard coded into the run_from_xml(). Leave out the febioCommand when running from the terminal.

LogPostProcessing.py
Python script creates a folder of images and csv files with model results. Can be run using Exp_to_Mod.xml with the run_all_in_file() function.

model_prediction_errors.py
To run after LogPostProcessing.py. Compares the kinematics between the model results and the experiment kinematics. Creates images, and files storing the rms error between them. Can be run using Exp_to_Mod.xml with the from_xml() function. see main function at end of script for examples of use.

transfer_med_groups.py
Transfers nodes and element sets from one mesh to another for a given component by treating the groups as a binary field and applying that binary field to the new mesh to find the correspondence of nodes and elements in the new mesh. This correspondence is used to select nodes and elements to generate groups with the new node and element sets.

StiffnessFromLog.py
Python script extracts rigid body reaction force displacement data from simulation output (.log).

ModelReduction_rigids3.py
Removes all contacts and ties for the specified component that needs to be removed and converts it into a rigid body.

edge_lengths.py
Python script to calculate average edge lengths, standard deviations, minimum and maximum edge lengths for a given volume mesh (.med).
Note: Full path of Geometry.feb file has to be provided in the FeBio_custom.feb file if using Febio 2.9 in Windows. If using Linux, there is no need to supply the whole path as long as the files are in the same folder.