Computational simulation of bone remodelling post reverse total shoulder arthroplasty

TL;DR

This paper presents a computational model simulating bone remodelling after reverse shoulder arthroplasty, helping understand how bone adapts to altered loading conditions to improve implant fixation.

Contribution

It introduces a novel finite element-based computational model that accounts for bone remodelling in response to shoulder implant procedures, validated with benchmark problems.

Findings

Model accurately predicts bone density changes post-surgery

Remodelling response varies with loading conditions

Provides insights for improving implant design

Abstract

Bone is a living material. It adapts, in an optimal sense, to loading by changing its density and trabeculae architecture - a process termed remodelling. Implanted orthopaedic devices can significantly alter the loading on the surrounding bone, which can have a detrimental impact on bone ingrowth that is critical to ensure secure implant fixation. In this contribution, a computational model that accounts for bone remodelling is developed to elucidate the response of bone following a reverse shoulder procedure for rotator cuff deficient patients. The physical process of remodelling is modelled using continuum scale, open system thermodynamics whereby the density of bone evolves isotropically in response to the loading it experiences. The fully-nonlinear continuum theory is solved approximately using the finite element method. The code developed to model the reverse shoulder procedure is validated using a series of benchmark problems.

Figures

Tables

Table 1. Table 1: Material properties of the scapula

Material Property	Value		Source
Young’s modulus	$E$	$=193{\text{ N mm}}^{-2}$	[5, 26]
Poisson’s ratio	$\nu$	$=0.25$	[5, 26]
Initial density	${\rho }_0^{\ast }$	$=0.6{\text{ g cm}}^{-3}$	[17]
Speed of adaptation	$c$	$=0.4{\text{ d mm}}^{-2}$	[33]
Attractor stimulus	${\Psi }_0^{\ast }$	$=0.01375{\text{ N mm}}^{-2}$	[33]
Mass conduction coefficient	$K_0$	$=0.05{\text{ d mm}}^{-2}$
Time step size	$\Delta t$	$=10$ d	[17]

Equations

$$D_t\widetilde{{\rho }_0}$$

$$\text{Div}\mathbf{P}$$

$$\widetilde{{\rho }_0}=\frac{{\rho }_0-{\rho }_0^{\ast }}{{\rho }_0^{\ast }}$$

$$\Psi (\widetilde{{\rho }_0},\mathbf{F})$$

$${\Psi }^{\text{neo}}$$

$$\mathbf{P}={\rho }_0^{\ast }\left[\widetilde{{\rho }_0}+1\right]\frac{\partial \Psi (\widetilde{{\rho }_0},\mathbf{F})}{\partial \mathbf{F}}.$$

$$\mathbf{R}=K_0\text{Grad}\widetilde{{\rho }_0},$$

$$R_0=\frac{c}{{\rho }_0^{\ast }}\left[{\left[\widetilde{{\rho }_0}+1\right]}^{-m}{\Psi }_0(\widetilde{{\rho }_0},\mathbf{F})-{\Psi }_0^{\ast }\right],$$

Full text

††journal: arXiv

1 Introduction

Reverse total shoulder arthroplasty is a surgical procedure to restore the integrity and functionality of the shoulder joint for rotator cuff deficient patients. The rotator cuff muscles stabilise the dynamic shoulder joint. In the reverse procedure, the humeral head is removed and a stem is placed in the central channel of the humerus with a polyethylene cup press-fitted into the top of the stem, see Figure 1. A hemispherical implant, termed the glenosphere, is screwed into the glenoid component of the scapula (shoulder blade). The deltoid muscle now compensates for the absent rotator cuff muscles, thereby restoring most of the functionality of the shoulder.

Even though the procedure restores functionality, many complications are reported post procedure - up to 75 % in some clinical series [3]. The most common complication is scapular notching [see 3, 4, 5, 6, 7], where, at full adduction, the humeral cup impinges on the inferior border of the scapula. The impingement causes bone resorption in the affected regions. Further, the polyethylene debris causes infection and osteolysis [4]. It is also hypothesised that the change in the stress distribution in the scapula post procedure may cause atrophy in the inferior region [8].

Computational models are a valuable, predictive tool in orthopaedics and for the design of implants. The effectiveness of different implant designs can be investigated initially with a computational model and, by following an iterative virtual design process, provide the patient with the optimal prosthesis.

The stress distribution in the bone affects the fixation quality of the implant. Bone is a complex, living tissue which responds to changes in the loading environment by altering its density and trabecular orientation, resulting in modified material properties. This phenomenon is termed bone remodelling. Bone consists of layers of dense cortical bone and spongy trabecular bone which is composed of strut-like structures surrounded by bone marrow. Furthermore, bone consists of a solid and a fluid phase. The remodelling response occurs in the solid phase and only at low strain rates. Remodelling is governed by cells in the solid phase, which deposit (osteoblastic action) and remove (osteoclastic action) material from the bone, resulting in density changes. For more details on the structure of bone and bone remodelling see [9, 10, 11].

Bone is often modelled as a porous solid. Because the remodelling response of bone is optimal to the loading environment, objective functions have been proposed to constrain and control the density evolution. The difference between a mechanical stimulus, which may be a stress, strain or strain-energy density, and an attractor state, which is the optimal state the bone is driving towards, provides the driving force for the density evolution, [see 11, 12, 13, for further details on objective function based bone remodelling models].

The open system thermodynamics approach introduced by Kuhl and Steinmann [14] is adopted here to describe bone remodelling in the scapula pre and post the reverse shoulder procedure. The authors are aware of only six works in the open literature on bone remodelling studies in the scapula [15, 16, 17, 18, 19, 20]. In these studies various remodelling theories were adopted that differ from the open system approach. Sharma and Robertson [17], for example, employ a strain-energy based, isotropic remodelling theory based on Jacobs et al. [12] together with a linear elastic material model. The initial geometry was provided by computed tomography (CT) scan data. The simulation results were reported to be over- and underestimating the actual density values determined from the CT scan. Similarly, Campoli et al. [19] used a CT scanned scapula as the model geometry. An optimisation based density evolution algorithm [21] with a strain-based mechanical stimulus was used. The muscle loads were obtained from the Delft Shoulder and Elbow model [22]. Campoli et al. compared the simulation results to the actual density values from the CT scan and reported an error less than 30%. Quental et al. [18] also used a three dimensional, subject-specific scapula as the geometry. Bone was modeled as a linear elastic, orthotropic, porous material following Fernandes et al. [23] [see also 24]. The optimisation algorithm involved a strain-based mechanical stimulus. The bone microstructure is captured by periodic unit cells accounting for the porosity [24], where the density is dependent on the solid volume fraction of the cells. The loads were obtained from a biomechanical model [18]. Quental et al. compared the density to the actual density values from the CT scan and reported an absolute error of below 33%.

The novel contribution of the work presented here is the application of the fully-nonlinear remodelling theory of Kuhl and Steinmann [14] to investigate density changes following reverse total shoulder arthroplasty. In addition, new insights into the theory are provided via a series of numerical examples. The validated finite element code developed to perform the numerical examples is available at [25].

This manuscript is organised as follows. The governing and constitutive equations for bone remodelling are detailed in Section 2. The numerical model is introduced in Section 3. The theory is then validated with a benchmark problem and different aspects of the model are investigated in Section 4. The validated code is applied to the intact scapula, followed by the application to the reverse shoulder in Section 5. Here, the ASTM F2028 [26] problem of a glenoid component implanted into a polyurethane foam block bone substitute is investigated. Finally, the reverse shoulder procedure is analysed. A discussion of the numerical results and model features are presented and conclusions are drawn in Section 6.

Notation

The scalar product of two second-order tensors $\bf{A}$ and $\bf{B}$ is denoted by $\bf{A}:\bf{B}$. The material time derivative of a quantity $\{\bullet\}$ is denoted by $D_{t}\{\bullet\}$. The divergence $\text{Div}\{\bullet\}$ and gradient $\text{Grad}\{\bullet\}$ operators are taken with respect to the material placement $\bf{X}$. As the formulations are based on the open system theory, a quantity $\{\bullet\}$ may be expressed per unit volume $\{\bullet\}_{0}$ or per unit mass $\{\bullet\}$, both with respect to the material configuration. The two descriptions are related via the material density $\rho_{0}$ as $\{\bullet\}_{0}=\rho_{0}\{\bullet\}$. $\bf{I}$ denotes the second-order identity tensor. For further details, see Kuhl and Steinmann [14].

2 Governing and constitutive equations for bone remodelling

Bone is modelled here as a continuum; i.e. is treated as continuous distribution of matter at the macroscopic scale. The microstructure of bone is thus not directly accounted for. A brief summary of the relevant nonlinear continuum mechanics is now given, for additional details see Holzapfel [27], for example.

Consider the deformation of a body $B_{0}\subset\mathbb{R}^{3}$ from the reference configuration to a current configuration occupied by the deformed body $B$, as shown in Figure 2. The control regions $\Omega_{0}\subset B_{0}$ and $\Omega\subset B$ with boundaries $\partial\Omega_{0}$ and $\partial\Omega$, respectively, are shown in Figure 2. The nonlinear deformation map $\boldsymbol{\varphi}$, i.e. the motion, maps the reference position $\bf{X}$ to the current position $\bf{x}$ at time $t$. The gradient of the nonlinear deformation map $\bf{F}:=\text{Grad}\boldsymbol{\varphi}$ is termed the deformation gradient. The determinant of the deformation gradient is defined by $J:=\text{det}\bf{F}$ and relates quantities such as the infinitesimal volume elements $dV$ and $dv$ between the two configurations. The right Cauchy–Green tensor ${\bf{C}}:={\bf{F}}^{T}{\bf{F}}$ is the chosen measure of deformation.

The theory developed here is based on open system thermodynamics, where the balance relations are formulated with material reference. For the sake of brevity, only the final governing equations and constitutive relations are stated. For a detailed derivation of the balance relations see Kuhl and Steinmann [14].

The balances of mass and linear momentum are given by

[TABLE]

The relative density $\tilde{\rho_{0}}$, introduced as the primary variable in the balance of mass (1), is defined by

[TABLE]

and $\rho_{0}^{*}$ is the initial density (assumed here to be uniform). $\bf{R}$ is the mass flux vector and $R_{0}$ is the scalar mass source term.

The balance of linear momentum (2) has been expressed in reduced form as the equilibrium equation by neglecting the inertial term. The time scales of the inertial and remodelling effects are significantly different, with remodelling taking place over a significantly longer time span. Body forces are also commonly neglected for clinical problems, as the applied loading during everyday activities is far greater than gravity and is therefore the stimulus for remodelling [12].

The constitutive relations proposed by Kuhl and Steinmann [14] are parameterised by the relative density $\tilde{\rho_{0}}$ and the right Cauchy–Green tensor $\bf{C}$. Bone is modelled as a compressible neo-Hookean material with a dependence on the density. The corresponding free energy is given by

[TABLE]

where $\lambda$ and $\mu$ are the Lamé constants, and $n$ is a material parameter which describes the bones porosity and is commonly chosen as 2 [28].

The constitutive relation for the first Piola–Kirchhoff stress tensor is given by

[TABLE]

The stress $\bf{P}$ is therefore governed by the coupled influence of the deformation and the density.

The mass flux vector $\bf{R}$ is defined analogous to Fick’s law of concentrations or Fourier’s law of heat conduction and is given by

[TABLE]

where $K_{0}\geq 0$ is the mass conduction coefficient.

The structure of the mass source term $R_{0}$ is based on the theory developed by Harrigan et al. [13]. Following Kuhl and Steinmann [14] the mass source term is given by

[TABLE]

where $\Psi_{0}=\rho_{0}\Psi$. $c$ governs the speed of the remodelling process and $m$ was introduced by Harrigan [29] to ensure numerical stability. $\Psi_{0}^{*}$ is the attractor stimulus, which is the reference free energy the system is driving towards.

3 Numerical model

The finite element method is used to approximately solve the nonlinear governing equations developed in the previous sections. The displacements and the relative density are introduced as the primary degrees of freedom. An implicit backward-Euler time stepping scheme is used for the time discretisation. The non-linear, coupled system is linearised following a Newton–Raphson approach and solved monolithically. The finite element code is developed using the AceGen library [30]. The full source code required to reproduce several of the benchmark problems presented in Section 4 is available at [25].

4 Benchmark problems

In order to validate the code and elucidate features of the formulation, a prototype example of a one-dimensional model problem, set out by Kuhl and Steinmann [14], is presented. A further new example investigates the influence of the time step size on the mass flux response.

4.1 An investigation of the mass source term

A specimen is loaded in tension by a time-dependent forcing function defined piecewise as shown in Figure 3(a). The specimen has unit length and a cross-sectional of $A=0.5\text{m}^{2}$, a Young’s modulus $E=1\text{N m}^{-2}$ and a Poisson’s ratio $\nu=0$. The mass flux $\bf{R}$ is set to zero to investigate the influence of the mass source term. The initial density $\rho_{0}^{*}=1\text{g cm}^{-3}$ and the attractor stimulus is set to $\Psi_{0}^{*}=1\text{N m}^{-2}$. The exponents $n$ and $m$ are set to $2$ and $3$, respectively. The time step size is chosen as $\Delta t=0.1\text{ s}$.

The evolution of the relative density and the displacement history of a point on the loaded face is shown in Figure 3 (b) and (c), respectively. The simulations match exactly the results of Kuhl and Steinmann [14]. The density and displacement change abruptly as the initial force and the subsequent stepwise increases of the force are applied and then converge to an equilibrium state.

The evolution of the volume specific, density-weighted free energy $\Psi_{0}$ and the biological stimulus $\left[\tilde{\rho_{0}}+1\right]^{-m}\Psi_{0}$ are shown in Figure 3 (d). The abrupt changes and the subsequent convergence towards an equilibrium state can be observed. The coupling between the density, displacement and the free energy is also evident. It is important to note that the biological stimulus will always drive towards the attractor stimulus $\Psi_{0}^{*}$, which equals 1 Pa in this example. The other quantities are increasing with the increasing force, as expected.

The influence of the attractor stimulus $\Psi_{0}^{*}$ on the density evolution is investigated further, as shown in Figure 4 (a). The initial decrease of density observed in 3 (b) is due to the high value of the attractor stimulus. There is an inverse relationship between the attractor stimulus and the density. This is explained by the nature of the mass source - as the attractor stimulus is smaller, the increase in density is naturally larger.

Finally, the influence of the time step size is investigated, as shown in Figure 4 (b). The time step size does not influence the density evolution significantly. The convergence rate is affected, but the steady state solution is essentially unchanged.

4.2 Investigation of the influence of the time step size and mass flux

The choice of time step size becomes important if a mass flux term is introduced. Consider now a unit size, three-dimensional specimen with a non-uniform initial density distribution, as shown in Figure 5. No mass source or mechanical loading is applied to the system. The density evolution is governed by the mass flux alone, i.e. (1) reduces to $D_{t}\tilde{\rho_{0}}=\text{Div}{\bf{R}}$. The mass flux term acts to smooth the solution, producing a uniform distribution of $\tilde{\rho_{0}}=0.5$ after 0.4 s if a time step size of $\Delta t=0.1$ s is used. However, for a smaller time step size of $\Delta t=0.01$ s, equilibrium is reached after 0.18 s. For a larger time step size of $\Delta t=1$ s, equilibrium is reached after 2 s. The same behaviour is observed for the problem of linear heat equation in a rigid conductor, where Fourier’s law governs the conduction. This is expected from the nature of the parabolic governing relations, but has not been commented on before in the context of bone remodelling.

5 Numerical prediction of remodelling in the scapula pre and post reverse shoulder procedure

The bone remodelling theory, validated in the previous section, is applied to the intact scapula and to the scapula post procedure. Prior to this, in order to elucidate features of the theory using a simpler geometry than the scapula, the implant is virtually implanted into a block. Finally, the full model of the scapula post procedure is developed.

5.1 Response of the scapula prior to implant

A geometry of the scapula produced by Mutsvangwa et al. [31] was used. The geometry is a statistical shape model, developed from 70 magnetic resonance imaging (MRI) scans. The pre-processing package ANSA [32] and the hexablocking technique are used to mesh the geometry using hexahedral elements.

The boundary conditions were approximated from those reported by Sharma and Robertson [17], who applied the muscle forces at 900 abduction, as shown in Figure 6. The levator scapula and rhomboid forces were omitted as these are relatively small. All forces were assumed to act normal to the surface.

The material properties for the scapula were chosen from various sources in the literature and are detailed in Table 1. The mass conduction coefficient is chosen as $K_{0}=0.05\text{ d mm}^{-2}$ to adjust the bounds of the density to fall within a physically reasonable range between $\rho_{0}=0\text{ gcm}^{-3}$ and the density of cortical bone $\rho_{0}=1.8\text{ gcm}^{-3}$ as reported by Sharma and Robertson [17].

The predicted density distribution is shown in Figure 7. The solution converges two time steps later than reported by Sharma and Robertson [17] at $t=$ 120 d. The resulting minimum and maximum relative density values of $\tilde{\rho_{0}}_{\text{min}}=$ -0.9734 and $\tilde{\rho_{0}}_{\text{max}}=$ 1.934 correspond to densities $\rho_{0}=0.016\text{ gcm}^{-3}$ and $\rho_{0}=1.76\text{ gcm}^{-3}$, respectively. The density distribution is similar to that reported by Sharma and Robertson [17]. The discrepancies are explained by the differences in geometry, forces, muscle attachment points, material parameters and, most importantly, the different bone remodelling theory.

5.2 Block example

This example investigates the remodelling response of a block that replicates the polyurethane foam block used in the ASTM F2028 [26] testing procedure to investigate glenoid component loosening. However, the block now has the potential to remodel like bone. A CAD model of the prosthesis was virtually implanted into the foam block. Owing to the symmetry of the problem, only half of the geometry was modelled and symmetry boundary conditions applied to the symmetry plane. The bottom edge was completely constrained and the left and right edges were constrained in the $x$-direction. A compressive load of 85 % of a standard body weight, which equals 750 N, was applied to the glenosphere. The mesh and boundary conditions are shown in Figure 8 (a).

The same material properties as prescribed for the scapula were used for the block (see Table 1). The glenosphere, screws and metaglene were modelled as linear elastic materials with a Young’s modulus of $E=230\text{ }000\text{ N mm}^{-2}$ for the glenosphere and $E=117\text{ }000\text{ N mm}^{-2}$ for the screws and metaglene. The Poisson’s ratio of the three linear elastic components was set to $\nu=0.3$. All components were modelled as perfectly bonded. The assumption that the density is continuous across the bone-screw interface leads to the artificial smoothing of the the density and related parameters as shown in Figure 8. This could be addressed by double defining the density degrees of freedom at the interface.

The density, and norms of the displacement and Cauchy stress $\boldsymbol{\sigma}:=J^{-1}\mathbf{P}\mathbf{F}^{T}$ distributions at $t=500$ d with $\Delta t=10$ d and $K_{0}=0.01\text{dmm}^{-2}$ are shown in Figure 8 (b), (c) and (d), respectively. As expected, the stress is concentrated in the screws as these are stiff, load bearing components. As a result, the block’s density increases in the vicinity of the screws. As the greater part of the block does not experience any loading, the density decreases in these areas. It is clear from the results that living bone would not evolve to this initial non-optimal cubic shape. Because the density is coupled to the material properties through the free energy, the bone is less stiff in the regions of low density and therefore offers little resistance to deformation.

The influence of the time step size and mass conduction coefficient on the convergence rate of the problem is now investigated. Increasing the conduction coefficient, tightens the bounds of the predicted density (i.e. the difference between the maximum and minimum density values throughout the domain) and also increases the rate of convergence to an equilibrium state. As shown for the benchmark problem in Figure 4 b, changing the time step size results in the same equilibrium state, however, the time to equilibrium is different. The time to equilibrium therefore has no physical interpretation.

5.3 Response of the scapula post reverse shoulder procedure

The scapula post reverse total shoulder arthroplasty is now investigated. The statistical shape model of the scapula is the one used in Section 5.1. The reverse shoulder prosthesis is virtually implanted in the scapula, using the same implant geometry as for the block example in the previous section. The surgical guidelines defined in the DePuy Delta XTEND reverse shoulder system’s surgical techniques [34] were followed and the virtual operation was performed in the software Mimics [35]. The various components of the model are shown in Figure 9 (a).

The magnitudes of the applied muscle forces were obtained from three patients whose shoulders were simulated using the Newcastle Shoulder Model [36] in OpenSim [37] - an open source musculoskeletal modeling software. The forces were also provided at 900 abduction, which is the same movement used for the intact scapula in Section 5.1. However, for the case of the reverse shoulder, the rotator cuff muscles (infraspinatus, supraspinatus, teres minor and subscapularis) are absent.

The Newcastle Shoulder Model outputs point-wise muscle attachment locations and forces for each muscle group. The point-wise forces for each muscle group are averaged per patient, and each of the muscle forces are averaged over the three patients. The muscle attachment regions were estimated from the attachment points in the Newcastle Shoulder Model, as well as from the BioDigitalHuman project [38]. The forces were again assumed to act normal to the surface. The muscle attachment regions and forces, as well as the fixed supports for the simulation are shown in Figure 9 (b).

The same material properties as for the intact scapula and the implant components in the block example were used. The time step size $\Delta t=10$ d and equilibrium was reached at $t=120$ d. The resulting density distributions are shown in Figure 10. The minimum and maximum relative density values of $\tilde{\rho_{0}}_{\text{min}}=-0.9745$ and $\tilde{\rho_{0}}_{\text{max}}=0.6651$ correspond to densities $\rho_{0}=0.0153\text{ g cm}^{-3}$ and $\rho_{0}=0.999\text{ g cm}^{-3}$, respectively. The maximum density is significantly less than for the intact scapula. Diminished muscle forces act on the scapula, and the glenohumeral joint contact force, which is the largest force applied on the system, is applied directly onto the prosthesis. This results in the stresses concentrating in the screws, as observed in the block example. As the screws are bearing the load, the stress in the bone is decreased, which results in a decrease in density. This is a known phenomenon in orthopaedics, termed stress shielding.

The density distribution in the vicinity of the implant is comparable to the results from the block example. In the regions where the muscle forces are applied, the density increases, as observed for the intact scapula. The density falls within the physiologically reasonable range.

It is of concern that the density decreases significantly following the reverse shoulder procedure in regions directly below and adjacent to the metaglene, as density loss here can lead to implant loosening and ultimately to the failure of the prosthesis, resulting in devastating consequences for the patient.

6 Discussion and conclusion

An open-system model for bone remodelling developed by Kuhl and Steinmann [14] was successfully applied to the case of a scapula post reverse shoulder arthroplasty.

The remodelling code developed was validated against a one-dimensional benchmark problem. Different aspects of the theory were investigated. It was shown that the attractor stimulus has to be small enough to ensure density apposition. The mass conduction coefficient $K_{0}$ should be calibrated to adjust the bounds of the density. The mass conduction coefficient and time step size were also shown to have an influence on the rate of convergence.

The validated model was applied to a statistical shape model of the scapula. The forces at 900 abduction were taken from literature and applied normal to the surface. It is difficult to compare the model to the literature as the geometry, muscle attachment points and even the bone remodelling theory are significantly different. Nevertheless, the results are reasonable once the mass conduction coefficient is calibrated so that physical density bounds are predicted. In most other remodelling theories, the density bounds are explicitly enforced, which is not necessary in the current model.

The shoulder was then examined post-procedure. First, the implant was placed in a block, which was modelled as bone, followed by an investigation of the scapula post procedure. The resulting density values in the scapula post procedure were significantly smaller than for the intact scapula, due to the absent rotator cuff muscles, and because the screws are the load bearing components in the system. The response to loading results in a densification around the screws and density resorption underneath and adjacent to the implant, which is a cause for concern. The areas where the muscle forces act, experience density apposition. A definite conclusion on the remodelling response can however not be drawn, as many uncertainties are present and only one loading scenario is considered.

In summary, the remodelling model provides physically reasonable and justifiable results. The density bounds are enforced naturally by adjusting the mass conduction coefficient to suit the problem. What is outstanding, is a comprehensive study of the scapula and calibration of the model parameters against experimental studies. Further, the exact geometry has to be imported to the muscoloskeletal shoulder model and a number of everyday loading scenarios should be considered. An extension of the model to include anisotropy considerations in the density evolution will be important as the density distribution in the scapula is anisotropic. Further, modelling of the contact between the implant and bone may yield a more predictive model. Once these extensions have been realised, the model provided here could be used as the basis for a reliable and powerful predictive tool in the research and design process of reverse shoulder arthroplasty.

Acknowledgements

The help and information for the reverse shoulder model provided by Jonathan Glenday from the Division of Biomedical Engineering at UCT and Andreas Kontaxis at the Leon Root Motion Analysis Laboratory at the Hospital for Special Surgery in New York City is greatly appreciated. Further, the statistical shape model provided by Tinashe Mutsvangwa in the Department of Biomedical Engineering at UCT is gratefully acknowledged.

H.L. and A.M. acknowledge funding from the National Research Foundation through the South African Research Chair in Computational Mechanics (SARCCM, grant agreement SARChI-CM- 47584) and the funding from the University of Cape Town through the Master’s Research Scholarship. The financial assistance of the National Research Foundation (NRF) through the Scarce Skills Master’s Scholarship (NRF Grant Number 102021) towards this research is hereby acknowledged. Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF.