Patient and imaging

We retrospectively recruited one participant (F, 19 years old) with unilateral LCPD deformity (Stulberg grade IV [2]), who had taken part in previous MRI studies [23, 24].

An MRI scan of both hip joints was obtained using a GE Discovery MR750 3 T scanner (Waukesha, WI, USA), with a MERGE water excitation sequence (field of view = 320 mm, slice thickness = 1.4 mm, slice spacing = 0.7 mm, acquisition matrix = 512 × 384, echo time = 16.1 ms, repetition time = 30.4 ms, flip angle = 5°).

Segmentation and mesh preparation

One rater (LJ) manually segmented the bone, articular cartilage, and labrum (including the transverse acetabular ligament) for each side using 3D Slicer v5.3.0 [25]. Segmentation labelmaps were then exported as surface mesh stereolithography files with 20 iterations of Laplacian “joint smoothing” to remove marching cubes artifacts while retaining watertight boundaries between adjacent surfaces.

Final smoothing and decimation of the chondrolabral meshes was performed in MeshLab v2025.07 [26], using the following steps:

1.

Quadric edge collapse decimation by 50%,

2.

Taubin smoothing,

3.

Isotropic explicit remeshing.

The collective aim of smoothing in MeshLab was to prepare the surface meshes for tetrahedralization by a) removing sharp edge features introduced by the joint smoothing algorithm, and b) encouraging faces with low aspect ratios (close to equilateral triangles). A full description of the parameters used for each step is available in Online Resource 1.

The target edge length for isotropic explicit remeshing was iteratively reduced by 10% until stress outputs changed by less than 5%, after which the mesh was considered converged. The final average edge length for the surface meshes was 0.72 mm.

ArtiSynth model

We performed simulation using the open-source Java-based ArtiSynth platform (v3.8) [27] and the project code is available in Online Resource 1. Bony structures were defined to be rigid bodies, and the articular cartilage and labrum were modelled using tetrahedral meshes generated within ArtiSynth from the smoothed surface meshes.

Boundary conditions

The pelvic bone was fixed, and the femur was free to move in all six degrees of freedom (Fig. 2). Nodes of FE meshes within 0.2 mm of the bone surface were considered part of the osteochondral boundary in order to accommodate vertex displacement during the mesh smoothing protocol. These nodes were attached to the bone, and frictionless contact was defined between the articulating femoral and acetabular finite element meshes (Fig. 2). Tests of typical friction coefficients for articular cartilage (µ = 0.001 and 0.01 [28, 29]) resulted in negligible changes to the model’s output at the expense of considerable increases in computation time (Online Resource 1).

Fig. 2

A simplified 2D diagram showing the boundary conditions and constraints in the digital twin model (affected side). The pelvic bone is fixed to the world coordinate frame. The femur is free to move with 6 degrees of freedom in response to force/torque inputs and constrained by frictionless contact (red) between the acetabular and femoral finite element meshes (blue) representing the joint’s articular cartilage and acetabular labrum

Joint rotations were input as anatomical rotations (flexion, adduction, and internal rotation) using a standard hip joint coordinate system [28] and we produced a class of methods to convert between anatomical angles and 3D rotation matrices in ArtiSynth. User-defined joint rotation trajectories were applied to the femur via a simple proportional-derivative torque control loop and the corresponding joint contact forces were applied directly to the femoral bone at the femoral head center (Fig. 2).

Material properties

We modeled articular cartilage as a Neo-Hookean hyperelastic material, with material parameters for dynamic and static loading derived from Park et al. (2004) and Athanasiou et al. (1995) respectively (Table 1) [20, 29].

Table 1 Neo-Hookean hyperelastic material parameters for articular cartilage and labrum matrix models

We represented the labrum as a transversely isotropic hyperelastic material region of the acetabular FE mesh, consisting of a Neo-Hookean matrix material base and a directional overlay representing the circumferential fiber orientation. We derived equilibrium matrix and fiber material parameters from Ferguson et al. (2001) [30]. In the absence of literature on the dynamic material properties of the acetabular labrum we assumed the matrix behaves similarly to articular cartilage, in that it has an instantaneous modulus 8 × greater than aggregate (equilibrium) modulus measured from the same sample (Table 1) [31]. We used custom Java code to calculate the fiber direction at each FE node from the labral geometry during the ArtiSynth model’s instantiation (Online Resource 1). The fiber overlay material was based on the passive properties of an ArtiSynth GenericMuscle material [32], with a toe region response described by an exponential stress coefficient of 0.051 MPa and uncrimping factor of 36, an exponential/linear transition stretch of 1.103, and a straightened fiber modulus of 74.7 MPa [30].

Simulation and outputs

Simulation was performed using ArtiSynth’s semi-implicit Newmark integrator with a simulation time step size of 0.002 s. FE meshes were represented with lumped-mass linear tetrahedral elements [27]. The output of the simulation was peak maximum-shear stress. At each sample time step (ts = 0.01 s), the stress for each FE node was calculated, and the highest values in the femoral cartilage and acetabular cartilage/labrum were recorded.

Dynamic gait simulation

We evaluated the dynamic performance of the digital twin model by simulating a typical level gait cycle in both hips, with the unaffected contralateral hip representing normal anatomy. We obtained joint contact force and angle trajectories for a typical level gait cycle from the literature [33, 34] and contact forces were scaled to the participant’s body weight (BW) for simulation.

Sensitivity analysis

We assessed the dynamic model’s sensitivity to expected variation in material properties and variation in joint angles due to motion capture error. We tested two material scenarios: one with a stiffer estimate of material parameters as used in some previous studies [14, 15] and one with more compliant properties representative of osteoarthritic cartilage [19] (Table 1). Additionally, we simulated four motion capture error scenarios, adding a constant ± 5° of flexion or ± 10° of internal rotation to the joint angle trajectory, representative of expected errors in marker-based motion capture systems during low-flexion activities such as gait [21, 22].

To test the repeatability of the baseline model, we ran it for two consecutive gait cycles and compared the peak stress output across the first and second cycles. To quantify differences between repeated gait cycles, between sensitivity analysis scenarios and between the affected and contralateral hips, we calculated the mean absolute percentage error (MAPE) for the whole gait cycle, separately for the femoral and acetabular meshes on the affected and contralateral sides:

$$\text = \frac\sum_^\frac_}-_}\right|}_}}$$

where \(n\) is the number of samples recorded. We also calculated the mean percentage error (MPE) to determine whether different scenarios resulted in a consistent over- or underestimate of peak stress compared to the baseline model:

$$\text = \frac\sum_^\frac_}-_}}_}}$$

The value of MPE always lies in the range ± MAPE. Positive and negative values MPE indicate test scenarios resulting in higher and lower peak stresses on average than the baseline simulation. If MPE =  + MAPE, the test scenario peak stress is higher than baseline at every sample time point, and lower at every time point if MPE = − MAPE.

Static high flexion simulation

We evaluated our digital twin model’s ability to accurately predict expected regions of high stress and femoral translations during impingement by reproducing static neutral and high flexion postures that had been imaged in an upright open MRI scanner (MROpen, ASG Superconductors, Genoa, Italy) [23].

The upright open MRI protocol included axial scans of the iliac crests and the distal femur in the supine posture, and sagittal and oblique scans of the hip joint in all postures. The axial and supine sagittal images provided the anatomical landmarks used to define each hip’s joint coordinate system [28]. We obtained transformation matrices for the pelvis and femur between each posture by rigid registration of the bones segmented from each sagittal image and converted these transformations to anatomical rotations for the simulation.

To simulate each posture (supine; seated; supine with high flexion, adduction, and internal rotation (FADIR); and seated FADIR), we first applied a constant force of 30% BW, typical of passive supine or seated postures [34], normal to the acetabular opening. The position and orientation of the oblique imaging planes, used to visualize anterior impingement in each posture in the upright open MRI, were obtained from image metadata. Custom MATLAB (version R2024a) code, available in Online Resource 1, was used to transform these planes to the starting femoral orientation, after which they were fixed to the femoral reference frame. The femur was then slowly rotated to the posture of interest and allowed to equilibrate, before taking a cross-section of the model in the oblique imaging plane for that posture and superimposing the cross-section stress map and bone outlines onto the corresponding MRI slice. From these composite images, we qualitatively assessed how well the model’s femoral and pelvic bone outlines matched the real-world MRI, and whether the model predicted stress concentrations corresponding with visible impingement locations.