Models

This module contains computational models relevant to bone remodelling. The Ruffoni Model contains the implementation of the bone mineralisation model by Ruffoni et al. (2008, 2009) to describe the evolution of BMDD curves over time depending on mineralization kinetics.

The Modiz Model contains the implementation of the bone mineralisation model by Modiz et al. (2026) to analyze porosity-dependent mineralization based on the boomerang-shaped curve, that is observed experimentally between apparent and material bone density.

Ruffoni Model

class bone_models.bone_mineralisation_models.models.ruffoni_model.Ruffoni_Model(simulation_time, number_of_grid_points, start='BMDD')[source]

Bases: object

Implements the Ruffoni et al. (2007) model for Bone Mineralization Density Distribution (BMDD) using FiPy’s Finite Volume Method for robust PDE solving.

Note

Source Publication: Ruffoni D., Fratzl P., Roschger P., Klaushofer K., Weinkamer R. (2007). The bone mineralization density distribution as a fingerprint of the mineralization process. Bone, 40, 1308–1319. DOI: 10.1016/j.bone.2007.01.012

The model solves the following partial differential equation:

\[\frac{\partial \rho(C, t)}{\partial t} + \frac{\partial \left( \rho(C, t) \, v_m(C) \right)}{\partial C} = - r(C, t) \, \rho(C, t)\]

with boundary condition:

\[\rho(0, t) \, v_m(0) = f(t)\]

Where:

  • \(\rho(C, t)\) — bone mineral density distribution (BMDD)

  • \(v_m(C)\) — mineralization velocity

  • \(r(C, t)\) — resorption rate at time \(t\)

  • \(f(t)\) — formation rate of new bone at time \(t\)

Initialize the Ruffoni BMDD model:

There are two options to initialize the steady state, determined by the start parameter:

  • ‘BMDD’ — Initialize from a reference BMDD distribution. The mineralization velocity and mineralization law are calculated from this steady-state BMDD using analytical formulas.

  • ‘mineralization law’ — Initialize from a predefined mineralization law. The BMDD and mineralization velocity are calculated from this mineralization law using analytical formulas.

Parameters:

  • simulation_time (float) – Total simulation time in years.

  • number_of_grid_points (int) – Number of grid points for the 1D calcium mesh.

  • start (str) – Initialization method, either 'BMDD' or 'mineralization law'.

Raises:

ValueError – If the start parameter is not 'BMDD' or 'mineralization law'.

During initialization, the model sets up a 1D mesh for calcium content, initializes the BMDD cell variable, and defines the mineralization law and velocity according to the specified starting condition. The mesh is defined over the range of calcium content from minimum_content to maximum_content as specified in the Ruffoni_Parameters class.

calculate_bone_volume()[source]

Calculate current total bone volume by integrating BMDD. The bone volume is calculated as the sum of the product of BMDD and cell volumes across all cells. The BMDD is averaged over the mesh cell volumes to ensure correct scaling.

Returns:

Total bone volume in micro m^3.

Return type:

float

calculate_formation_rate(t)[source]

Return the formation rate at time t.

Parameters:

t (float) – Time in years.

Returns:

Formation rate at time t.

Return type:

float

calculate_resorption_rate(t)[source]

Return resorption rate at time t.

Parameters:

t (float) – Time in years.

Returns:

Resorption rate at time t.

Return type:

float

check_mass_conservation(formation_rate, resorption_rate, bone_volume)[source]

This function checks the mass conservation in PDE solving - should be close to zero for equilibrium.

Parameters:

  • formation_rate (float) – Formation rate at current time step.

  • resorption_rate (float) – Resorption rate at current time step.

  • bone_volume (float) – Current total bone volume.

Returns:

Change in bone volume (should be close to zero for mass conservation).

Return type:

float

initialize_BMDD_from_mineralization_law(plot=False)[source]

Initialize the Bone Mineralization Density Distribution (BMDD) from a prescribed mineralization law. This method is called when the model is initialized with start='mineralization law'.

The initialization is based on the steady-state analytical formulation described in the source paper. The BMDD is computed by integrating the ratio of the resorption rate to the mineralization velocity across the calcium content range, taking the exponential of the integral and multiplying by a factor depending on formed volume and velocity. The resulting BMDD is then smoothed using a Gaussian filter and interpolated to match the mesh cell centers used in the numerical model.

If plot=True, the function additionally displays: - The exponential of the integrated inverse mineralization law (exp(-∫ r/v)), - The initial smoothed BMDD distribution over calcium content.

Parameters:

plot (bool) – If True, plots the exponential of the inverse mineralization integral and the smoothed initial BMDD distribution. Defaults to False.

Returns:

The initialized BMDD distribution mapped to the model mesh.

Return type:

numpy.ndarray

Details: - The BMDD is initialized as: \(\rho(C) = \frac{f(0)}{v_m(C)} \exp\left(-\int_{C_{min}}^{C} \frac{r(C')}{v_m(C')} \, dC'\right)\) where:

  • \(f(0)\) is the formed bone volume at steady state,

  • \(r(C)\) is the resorption rate,

  • \(v_m(C)\) is the mineralization velocity.

  • The integral is evaluated numerically using adaptive quadrature.

  • The smoothed BMDD is interpolated with a monotone cubic spline (PchipInterpolator) and clipped to enforce non-negativity.

initialize_BMDD_from_reference(plot=False)[source]

Initialize the Bone Mineralization Density Distribution (BMDD) from a reference skewed Gaussian distribution. This method is called when the model is initialized with start='BMDD'.

The reference BMDD parameters are taken from Roschger et al. (2003), providing a physiological steady-state distribution representative of healthy trabecular bone. The distribution is constructed using a skew-normal probability density function (PDF) defined by its peak calcium content and spread, and is scaled to match the reference peak amplitude (×25 for unit consistency).

If plot=True, the method plots the initial BMDD distribution over calcium content.

Parameters:

plot (bool) – If True, plots the initial BMDD distribution. Defaults to False.

Returns:

The initialized BMDD distribution defined over the mesh calcium content range.

Return type:

numpy.ndarray

Details: - The reference BMDD is computed as: \(\rho(C) = A \cdot \text{SkewNorm}(C; a, \mu, \sigma)\) where:

  • \(a\) is the skewness parameter (negative for left-skew, positive for right-skew),

  • \(\mu\) is the peak calcium content (from Roschger reference data),

  • \(\sigma\) is the mean standard deviation of lower and upper percentiles.

  • The resulting distribution is normalized such that its integral equals 1 and scaled by a constant factor (25) to match the reference BMDD peak.

  • The function optionally visualizes the initialized BMDD.

initialize_double_exponential_mineralization_law(time_values)[source]

Initialize a double-exponential mineralization law, which maps tissue age (time since formation) to calcium content. This function evaluates the mineralization law at the provided time_values.

The mineralization law is defined as:

\[C(t) = C_1 (1 - \exp(\frac{t}{r_1}))+ (C_2 - C_1) (1 - \exp(\frac{t}{r_2}))\]

where:

  • \(C_1, C_2\) are the asymptotic calcium contents of the fast and slow mineralization phases,

  • \(r_1, r_2\) are characteristic time constants of the two phases,

  • \(t\) is the time/ tissue age.

Parameters:

time_values (np.ndarray) – Time array at which the mineralization law is evaluated.

Returns:

Calcium content values corresponding to the provided time points.

Return type:

np.ndarray

initialize_double_hyperbolic_mineralization_law(time_values)[source]

Initialize a double-hyperbolic mineralization law, which maps tissue age (time since formation) to calcium content. This function evaluates the mineralization law at the provided time_values.

The mineralization law is defined as:

\[C(t) = C_1 \frac{t / r_1}{1 + (t / r_1)} + C_2 \frac{t / r_2}{1 + (t / r_2)}\]

where:

  • \(C_1, C_2\) are the asymptotic calcium contents of the fast and slow mineralization phases,

  • \(r_1, r_2\) are characteristic time constants of the two phases,

  • \(t\) is the time/ tissue age.

Parameters:

time_values (np.ndarray) – Time array at which the mineralization law is evaluated.

Returns:

Calcium content values corresponding to the provided time points.

Return type:

np.ndarray

initialize_formed_bone_volume(bone_volume, start)[source]

Initialize the formed bone volume based on the initial BMDD and initial bone volume. This is used to set the initial condition for the BMDD.

Parameters:

  • bone_volume (float) – The total bone volume calculated from the initial BMDD.

  • start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law

Raises:

ValueError – If start parameter is not ‘BMDD’ or ‘mineralization law’.

Returns:

None

initialize_inverse_mineralization_law(start, plot=False)[source]

Initialize the inverse mineralization law based on either BMDD or mineralization law. This method sets up the inverse mineralization law as a PchipInterpolator based on the calcium content and corresponding inverse mineralization law values.

Parameters:

  • start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law’.

  • plot (bool) – If True, plots the inverse mineralization law.

Raises:

ValueError – If start parameter is not ‘BMDD’ or ‘mineralization law’.

Returns:

None

initialize_inverse_mineralization_law_from_BMDD()[source]

Initialize the inverse mineralization law from BMDD steady state by integrating 1 / mineralization velocity from 0 to the calcium content value using numerical integration.

Returns:

Tuple of numpy arrays containing calcium values and corresponding inverse mineralization law values.

Return type:

tuple of np.ndarray

Raises:

ValueError – If the initial bone volume is less than or equal to zero.

initialize_inverse_mineralization_law_from_mineralization_law()[source]

Initialize the inverse mineralization law from the mineralization law by generating a range of time values and calculating the corresponding calcium values using the mineralization law. The method ensures that the calcium values are unique and monotonic, and clamps them to the valid range defined by the minimum and maximum calcium content.

Returns:

Tuple of numpy arrays containing calcium values and corresponding inverse mineralization law values.

Return type:

tuple of np.ndarray

initialize_mineralization_law(start, plot=False)[source]

Initialize the mineralization law based on either BMDD or given mineralization law. For the ‘BMDD’ start, it inverts the values of the inverse mineralization law. For ‘mineralization law’ start, a double hyperbolic or double exponential mineralization law is prescribed. Duplicates are removed, and the values are ensured to be monotonic. The values are interpolated using PchipInterpolator for smoothness.

Parameters:

  • start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law’.

  • plot (bool) – If True, plots the mineralization law. Defaults to False.

Raises:

ValueError – If start parameter is not ‘BMDD’ or ‘mineralization law’.

Returns:

None

initialize_mineralization_velocity(start, plot=False)[source]

Initialize mineralization velocity based on the initial BMDD or mineralization law by calling the respective functions. This method sets up the mineralization velocity law as a PchipInterpolator based on the calcium content and corresponding mineralization velocity values.

Parameters:

  • start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law’.

  • plot (bool) – If True, plots the mineralization velocity law.

Raises:

ValueError – If start parameter is not ‘BMDD’ or ‘mineralization law’.

Returns:

None

initialize_mineralization_velocity_from_BMDD()[source]

Initialize mineralization velocity from BMDD steady state when start is ‘BMDD’. The velocity is calculated by integrating the product of BMDD and resorption rate over the calcium content range and dividing the result by the corresponding BMDD value.

The integration is performed using the PchipInterpolator for smooth interpolation of BMDD and resorption rates and numpy quad for numerical integration.

Returns:

Mineralization velocity values as a numpy array and corresponding calcium values as a numpy array.

Return type:

tuple of np.ndarray

initialize_mineralization_velocity_from_mineralization_law(start, plot=False)[source]

Initialize mineralization velocity from the mineralization law when start is ‘mineralization law’. The mineralization velocity is calculated as 1 / derivative of the inverse mineralization law. This function calls the initialize_inverse_mineralization_law method to ensure the inverse law is set up correctly. The calcium values are generated over a range defined by the minimum and maximum calcium content. The PChipInterpolator.derivative() method is used to calculate the derivative. Extrapolation handled by clamping the values to avoid issues with division by zero or negative derivatives.

Parameters:

  • start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law’.

  • plot (bool) – If True, plots the derivative of the inverse mineralization law.

Returns:

Mineralization velocity values as a numpy array and corresponding calcium values as a numpy array.

Return type:

tuple of np.ndarray

initialize_model()[source]

Initialize the model based on the specified starting condition. If start is ‘BMDD’, it initializes a reference BMDD - in this case a Gaussian distribution - calculates the mineralization velocity, inverse mineralization law, and finally the mineralization law. If start is ‘mineralization law’, it initializes the mineralization law first, calculates the mineralization velocity, and finally the steady-state BMDD.

Raises:

ValueError – If the start parameter is not ‘BMDD’ or ‘mineralization law’.

initialize_steady_state_BMDD(start)[source]

Initialize BMDD depending on the start parameter.

If start is ‘BMDD’, it initializes from a reference BMDD distribution.

If start is ‘mineralization law’, it initializes from a mineralization law.

The BMDD is set to the initial value, and the formation rates are calculated based on the initial BMDD (formation depends on bone volume).

Parameters:

start (str) – Initialization method, either ‘BMDD’ or ‘mineralization law’.

Raises:

ValueError – If start parameter is not ‘BMDD’ or ‘mineralization law’.

Returns:

None

plot_results(BMDD_evolution, BV_evolution, time_points)[source]

Plot normalized BMDD and bone volume evolution at saved / specified time points. It also plots the mineralization law, mineralization velocity, and the formed and resorbed bone volume over time.

Parameters:

  • BMDD_evolution (np.ndarray) – Evolution of BMDD over time as a numpy array.

  • BV_evolution (np.ndarray) – Evolution of bone volume over time as a numpy array.

  • time_points (np.ndarray) – Time points at which the BMDD and bone volume were saved.

Returns:

None

solve_for_BMDD(save_interval=0.1)[source]

Solve the time evolution of the Bone Mineralization Density Distribution (BMDD) using FiPy’s Finite Volume Method (FVM).

The model is solved iteratively in time, updating velocity, resorption, and formation at each step. A transient term models time evolution, a convection term models advection by the mineralization velocity, and an implicit source term models resorption.

Boundary conditions: - Left boundary (low calcium): formation flux, given by the formed bone volume divided by the mineralization velocity at zero calcium. - Right boundary (high calcium): absorbing boundary (no-flux condition).

The PDE is solved using FiPy’s sweep method with adaptive time-stepping controlled by a residual constraint. Results are stored at regular intervals.

Parameters:

save_interval (float, optional) – Time interval (in years) between stored simulation states. Defaults to 0.1.

Returns:

(BMDD_evolution, BV_evolution, time_points) as BMDD, bone volume, and time points for each saved time step.

Return type:

tuple of numpy.ndarray

Raises:

ValueError – If the adaptive time step becomes too small during the simulation.

Notes

  • The velocity field is defined as a face variable because it represents BMDD flux across cell boundaries.

  • The resorption coefficient is a cell variable since it is a local property.

  • The adaptive solver halves the time step until the desired residual is achieved, or terminates if the minimum time step is reached.

Modiz Model

class bone_models.bone_mineralisation_models.models.modiz_model.Modiz_Model(load_case, porosity, specific_surface_multiplier=1, hypothesis=3)[source]

Bases: Martinez_Reina_Model

Implementation of the Modiz model for explicit bone mineralisation.

Based on the work of Modiz et al. (2026), this model simulates the mineralisation process in bone tissue as a function of porosity. It integrates bone ageing, volume fraction dynamics, turnover rates, and specific surface area to track mineral density over time.

Note

Source Publication: Corinna Modiz, Natalia M. Castoldi, Jose L. Calvo-Gallego, Stefan Scheiner, Vittorio Sansalone, Saulo Martelli, Javier Martínez-Reina, Peter Pivonka, Factors of mineralization dynamics — A mathematical model for microstructural bone remodeling Bone, 209, 117905. DOI: 10.1016/j.bone.2026.117905

Note

This model is a specialized subclass of Martinez-Reina Model and requires a Lerebours Model V1 (Legacy) instance to handle cell population steady-state calculations.

Key Components:

  • Ageing Queue: Tracks the temporal evolution of bone patches.

  • Porosity Dependence: Mineralisation kinetics scale with bone volume fraction.

  • Hypothesis Framework: Supports multiple mechanistic assumptions (1, 2.1, 2.2, 3) regarding resorption bias and mineral apposition.

Initialise the Modiz model with biological and geometric parameters.

Parameters:

  • load_case (Lerebours_Load_Case) – Load case defining steady-state cell concentrations.

  • porosity (float) – Initial porosity of the bone (percentage, e.g., 20.0).

  • specific_surface_multiplier (float) – Scaling factor for the specific surface area, used for dispersion/boundary curves.

  • hypothesis (float) – The model hypothesis to apply (1 (turnover), 2.1 (age bias), 2.2 (reduced resorption), or 3 (mineralization kinetics)).

calculate_average_mineral_content(OBa, OCa, bone_volume_fraction)[source]

Calculate the average mineral content based on the ageing queue after initialization, current bone volume fraction and mineralization law.

Parameters:

  • OBa (float) – active osteoblast cell concentration

  • OCa (float) – active osteoclast cell concentration

  • bone_volume_fraction (float) – bone volume fraction

Returns:

average mineral content

Return type:

float

calculate_error_for_turnover(fitting_parameter, bone_volume_fraction_data, experimental_turnover_data)[source]

Calculate the predicted turnover and then the error to the experimental turnover data for the specific bone volume fraction. The error is calculated as the sum of squared differences between the predicted and experimental turnover data.

Parameters:

  • fitting_parameter (float) – scaling parameter for the specific surface

  • bone_volume_fraction_data (numpy array) – bone volume fraction data

  • experimental_turnover_data (numpy array) – experimental turnover data

Returns:

error

Return type:

float

calculate_mineral_apposition_rate(bone_volume_fraction)[source]

Calculate the mineral apposition rate based on the initial bone volume fraction and prescribed hypothesis.

Parameters:

bone_volume_fraction (float) – initial bone volume fraction

Returns:

mineral apposition rate

Return type:

float

calculate_mineral_densities(bone_volume_fraction)[source]

Calculate the mineral densities based on the initial bone volume fraction and prescribed hypothesis. This function calculates the steady-state concentrations based on the Lerebours model, initialises the ageing queue and calculates the material and apparent density. Note: we use the turnover curve the other way around compared to the Lerebours model, so we need to put in bone volume fraction instead of porosity.

Parameters:

bone_volume_fraction (float) – initial bone volume fraction

Returns:

bone material density, bone apparent density

Return type:

tuple of floats

calculate_mineral_queue(ageing_queue, bone_volume_fraction)[source]

Calculate the mineral queue based on the ageing queue after initialization and initial bone volume fraction. The queue is calculated by multiplying the ageing queue with mineral content from the mineralisation law.

Parameters:

  • ageing_queue (numpy array) – ageing queue

  • bone_volume_fraction (float) – initial bone volume fraction

Returns:

mineral queue

Return type:

numpy array

calculate_mineralisation_law(t, bone_volume_fraction)[source]

Calculate the mineralisation law based on the age of the bone patch (queue element) and initial bone volume fraction. The law is used to calculate the mineral content in the bone tissue and includes lag time and exponential increase to the maximum mineral content. The mineral apposition rate depends on the initial bone volume fraction and prescribed hypothesis.

Parameters:

  • t (int) – age of the bone patch (queue element)

  • bone_volume_fraction (float) – initial bone volume fraction

Returns:

mineral content

Return type:

float

fit_specific_surface_to_turnover_data(show_plot=False)[source]

Fit the specific surface to the turnover data from Parfitt et al. The specific surface is calculated using the Lerebours model and the turnover data is fitted using a least-squares optimization and error function.

Parameters:

show_plot (bool) – whether to show and save the plot of the fitted curve, data and specific surface

Returns:

scaling parameter

Return type:

float

initialise_ageing_queue(OBa, OCa, bone_volume_fraction)[source]

Initialise the ageing queue with the initial bone volume fraction and resorbed/ formed fractions. It is initialised with zeros and the update_ageing_queue function is called until the queue is filled.

Parameters:

  • OBa (float) – active osteoblast cell concentration

  • OCa (float) – active osteoclast cell concentration

  • bone_volume_fraction (float) – bone volume fraction

Returns:

None

plot_turnover_fit(fitted_scaling_parameter, parfitt_bone_volume_fraction, parfitt_turnover)[source]

Plot the fitted turnover curve along with the experimental data and specific surface.

Parameters:

  • fitted_scaling_parameter (float) – fitted scaling parameter for the specific surface

  • parfitt_bone_volume_fraction (numpy array) – bone volume fraction data from Parfitt et al.

  • parfitt_turnover (numpy array) – turnover data from Parfitt et al.

Returns:

None

update_ageing_queue(OBa, OCa, bone_volume_fraction)[source]

Update the bone ageing queue for a single time step. This method shifts the queue applies resorption and formation rates. The update follows these steps:

  1. Shift & Formation: All bone patches age by one day; new bone (formed fraction) is added to the start of the queue.

  2. Resorption: Existing bone volume is reduced across the queue based on the active osteoclast (OCa) concentration.

  3. Hypothesis Handling:

    • If hypothesis 2.1: Resorption rate is biased by age.

    • If hypothesis 2.2: Resorption rate is scaled with the bone volume fraction.

    • Otherwise: Resorption follows standard lag-time rules.

  4. Mass Balance: The final queue element acts as a residual to ensure the total sum matches the current bone_volume_fraction.

Parameters:

  • OBa (float) – Active osteoblast cell concentration.

  • OCa (float) – Active osteoclast cell concentration.

  • bone_volume_fraction (float) – Current fraction of bone volume (0.0 to 1.0).

Notes:

This method modifies self.ageing_queue in-place.