Models

This module contains computational models relevant to bone remodelling. The Base Model is the base class for all models and contains the most general representation of a bone cell population models.

The Lemaire Model contains the bone cell population model based on the work of Lemaire et al. (2004) and is a submodule of the Base Model.

The Pivonka Model contains the bone cell population model based on the work of Pivonka et al. (2008) and is a submodule of the Lemaire Model. It adds more details to the intracellular communication pathways.

The Scheiner Model contains the bone cell population model based on the work of Scheiner et al. (2013) and is a submodule of the Pivonka Model. It adds the effect of mechanical loading on the bone cell population model using a homogenisation framework to calculate microscopic strain energy density from macroscopic stress tensors.

The Modiz Model contains the bone cell population model based on the work of Modiz et al. (2025) and is a submodule of the Lemaire Model. It changes the activation of osteoblasts by parathyroid hormone (PTH) to include details about pulse characteristics and two states of the receptor. The extension is based on the Martonova Model.

The Martonova Model calculates cellular activity constants using pulsatile PTH based on the work of Martonova et al. (2023).

The Martinez-Reina Model contains the bone cell population model based on the work of Martinez-Reina et al. (2019) and is a submodule of the Scheiner Model. It includes a bone cell population model for trabecular bone with mechanical loading, mineralisation, postmenopausal osteoporosis, and denosumab treatment.

The Lerebours Model contains the bone cell population model based on the work of Lerebours et al. (2015) and is a submodule of the Scheiner Model. It includes a bone cell population model with cell concentrations depending on the bone volume fraction and surface availability. It also includes a mechanical feedback mechanism based on microscopic strain energy density and bone site.

Base Model

class bone_models.bone_cell_population_models.models.base_model.Base_Model

Bases: object

This class implements the base model for bone cell population models. It essentially includes the ODE system that is shared across future models and the functions to solve for steady-state and solution. It is not intended to be used as a model itself but only indicates the most basic version of the bone cell population models. Functions that are specific to a certain model should be implemented in the respective model class, otherwise the functions return 0 (additive) or 1 (multiplicative) values.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• OCa (float) – osteoclast cell concentration

• parameters (Parameters) – Model parameters, see Parameters for details.

• initial_guess_root (numpy.ndarray) – initial guess for the root-finding algorithm for steady-state

• steady_state (object) – steady state values of the model

• t (float) – time variable

• tspan (numpy.ndarray with start and end time) – time span for the ODE solver

• x (list) – state variables of the model

• solution (scipy.integrate._ivp.ivp.OdeResult) – solution of the ODE system

• initial_bone_volume_fraction (float) – initial bone volume fraction

• dOBpdt (float) – rate of change of precursor osteoblast cell concentration

• dOBadt (float) – rate of change of active osteoblast cell concentration

• dOCadt (float) – rate of change of osteoclast cell concentration

• dxdt (list) – rate of change of state variables

• bone_volume_fraction (list) – bone volume fraction over time

Constructor for the Base_Model class.

apply_mechanical_effects(OBp, OBa, OCa, t)

Apply mechanical effects to the bone cell population model. Returns 0 (additive) as neutral value if not relevant to the specific model.

Returns:

mechanical effects acting on the bone cell population model

Return type:

float

apply_medication_effects_OBa()

Apply medication effects to active osteoblasts. Returns 1 (multiplicative) as neutral value if not relevant to the specific model.

Returns:

medication effects acting on active osteoblasts

Return type:

float

apply_medication_effects_OBp()

Apply medication effects to precursor osteoblasts. Returns 0 (additive) as neutral value if not relevant to the specific model.

Returns:

medication effects acting on precursor osteoblasts

Return type:

float

bone_cell_population_model(x, t=None)

Calculates the system of ordinary differential equations for the bone cell population model. This function is inherited by the specific models. If a function is not relevant in the specific model (e.g. mechanical effects), it returns 0 (additive) or 1 (multiplicative) as neutral values.

Parameters:

• x (list) – state variables of the model

• t (float) – time variable

Returns:

rate of change of state variables

Return type:

list

calculate_RANKL_activation_OCp(OBp, OBa, t)

Calculate the activation of precursor osteoclasts by RANKL.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

activation of precursor osteoclasts by RANKL

Return type:

float

calculate_TGFb_activation_OBu(OCa, t)

Calculate the activation of uncommitted osteoblasts by TGFb.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

activation of uncommitted osteoblasts by TGFb

Return type:

float

calculate_TGFb_activation_OCa(OCa, t)

Calculate the activation of active osteoclasts by TGFb.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

activation of active osteoclasts by TGFb

Return type:

float

calculate_TGFb_repression_OBp(OCa, t)

Calculate the repression of precursor osteoblasts by TGFb.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

repression of precursor osteoblasts by TGFb

calculate_bone_volume_fraction_change(solution, steady_state, initial_bone_volume_fraction)

Calculate the change in bone volume fraction over time.

Parameters:

• solution (scipy.integrate._ivp.ivp.OdeResult) – solution of the ODE system

• steady_state (numpy.ndarray) – steady state values of the model

• initial_bone_volume_fraction (float) – initial bone volume fraction

Returns:

bone volume fraction over time

Return type:

list

calculate_external_injection_OBa(t)

Calculate the external injection of active osteoblasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of active osteoblasts

Return type:

float

calculate_external_injection_OBp(t)

Calculate the external injection of precursor osteoblasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of precursor osteoblasts

Return type:

float

calculate_external_injection_OCa(t)

Calculate the external injection of active osteoclasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of active osteoclasts

Return type:

float

calculate_steady_state()

Calculate the steady state of the bone cell population model using root finding of the ODE system.

Returns:

steady state values of the model

Return type:

numpy.ndarray

solve_bone_cell_population_model(tspan)

Solve the bone cell population model using the ODE system over a given time interval. The initial conditions are set to the steady-state values.

Parameters:

tspan (numpy.ndarray with start and end time) – time span for the ODE solver

Returns:

solution of the ODE system

Return type:

scipy.integrate._ivp.ivp.OdeResult

Lemaire Model

class bone_models.bone_cell_population_models.models.lemaire_model.Lemaire_Model(load_case)

Bases: Base_Model

This class implements the bone cell population model by Lemaire et al. (2004) as a subclass of the Base_Model class.

Note

Source Publication: Lemaire, V., Tobin, F. L., Greller, L. D., Cho, C. R., & Suva, L. J. (2004). Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. Journal of Theoretical Biology, 229(3), 293-309. DOI: 10.1016/j.jtbi.2004.03.023

Parameters:

• load_case (object) – load case for the model

• parameters (Lemaire_Parameters) – model parameters

• initial_guess_root (numpy.ndarray) – initial guess for the root-finding algorithm for steady-state

• steady_state (object) – steady state values of the model

• t (float) – time variable

• tspan (numpy.ndarray with start and end time) – time span for the ODE solver

• x (list) – state variables of the model

• solution (scipy.integrate._ivp.ivp.OdeResult) – solution of the ODE system

• initial_bone_volume_fraction (float) – initial bone volume fraction

• dOBpdt (float) – rate of change of precursor osteoblast cell concentration

• dOBadt (float) – rate of change of active osteoblast cell concentration

• dOCadt (float) – rate of change of osteoclast cell concentration

• dxdt (list) – rate of change of state variables

• bone_volume_fraction (list) – bone volume fraction over time

Constructor for the Lemaire_Model class.

Parameters:

load_case (object) – load case for the model

calculate_OPG_concentration(OBp, t)

Calculate the OPG concentration depending on intrinsic and external OPG and degradation rate.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• t (float) – time variable

Returns:

OPG concentration

Return type:

float

calculate_PTH_activation_OB(t)

Calculate the activation of osteoblasts by PTH.

Parameters:

t (float) – time variable

Returns:

activation of osteoblasts by PTH

Return type:

float

calculate_PTH_concentration(t)

Calculate the PTH concentration depending on intrinsic and external PTH and degradation rate.

Parameters:

t (float) – time variable

Returns:

PTH concentration

Return type:

float

calculate_RANKL_activation_OCp(OBp, OBa, t)

Calculate the activation of precursor osteoclasts by RANKL.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

activation of precursor osteoclasts by RANKL

Return type:

float

calculate_TGFb_activation_OBu(OCa, t)

Calculate the activation of uncommitted osteoblasts by TGF-beta.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

activation of uncommitted osteoblasts by TGF-beta

Return type:

float

calculate_TGFb_activation_OCa(OCa, t)

Calculate the activation of active osteoclasts by TGF-beta.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

activation of active osteoclasts by TGF-beta

Return type:

float

calculate_TGFb_repression_OBp(OCa, t)

Calculate the repression of precursor osteoblasts by TGF-beta.

Parameters:

• OCa (float) – osteoclast cell concentration

• t (float) – time variable

Returns:

repression of precursor osteoblasts by TGF-beta

Return type:

float

calculate_bone_volume_fraction_change(time, solution, steady_state, initial_bone_volume_fraction)

Calculate the bone volume fraction over time numerically using the solution of the ODE system and the trapezoidal rule.

Parameters:

• time (numpy.ndarray) – time variable

• solution (numpy.ndarray) – solution of the ODE system

• steady_state (numpy.ndarray) – steady state values of the model

• initial_bone_volume_fraction (float) – initial bone volume fraction

Returns:

bone volume fraction over time

Return type:

list

calculate_external_injection_OBa(t)

Calculate the external injection of active osteoblasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of active osteoblasts

Return type:

float

calculate_external_injection_OBp(t)

Calculate the external injection of precursor osteoblasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of precursor osteoblasts

Return type:

float

calculate_external_injection_OCa(t)

Calculate the external injection of active osteoclasts (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of active osteoclasts

Return type:

float

calculate_external_injection_OPG(t)

Calculate the external injection of OPG (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of OPG

Return type:

float

calculate_external_injection_PTH(t)

Calculate the external injection of PTH (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of PTH

Return type:

float

calculate_external_injection_RANKL(t)

Calculate the external injection of RANKL (used in load case scenarios).

Parameters:

t (float) – time variable

Returns:

external injection of RANKL

Return type:

float

Martonova Model

class bone_models.bone_cell_population_models.models.martonova_model.Martonova_Model(load_case)

Bases: object

This class defines the two-state receptor model by Martonova et al. (2023) for pulsatile endogenous PTH in healthy and disease state including drug administration.

Note

Source Publication: Martonova D., Lavaill M., Forwood M.R., Robling A., Cooper D.M.L., Leyendecker S., Pivonka P. (2023). Effects of PTH glandular and external dosing patterns on bone cell activity using a two-state receptor model —Implications for bone disease progression and treatment PLOS ONE, 18(3), e0283544. DOI: 10.1371/journal.pone.0283544

It is based on the model by Li and Goldbeter (1986) for receptor-ligand binding, which will be referred to in some of the implemented model functions.

Li Y., Goldbeter A. (1989). Frequency specificity in intercellular communication. Influence of patterns of periodic signaling on target cell responsiveness. Biophysical journal, 55(1), 125–145. DOI: 10.1016/S0006-3495(89)82785-7

Parameters:

• load_case (Martonova_Load_Case) – Load case for the model

• parameters (Martonova_Parameters) – Model parameters

• initial_condition (numpy.ndarray) – Initial condition for the ODE system for receptor-ligand binding.

• number_of_periods (int) – Number of periods for basal PTH pulse.

• period_for_activity_constants (int) – Period for calculating activity constants after cellular adaptation to basal PTH.

Constructor method for the Martonova_Model class. It calculates the active complex activity constant and the drug PTH pulse based on the load case.

Parameters:

load_case (object) – Load case for the model, see Martonova_Load_Case for details.

PK_PD_model(t, x)

This function defines the pharmacokinetics-pharmacodynamics model based on an ODE system for drug PTH pulse. The ODE system describes the change of dose, drug concentration, and area under the curve over time depending on absorption, elimination, distribution rates and bioavailability.

Parameters:

• t (float) – time variable

• x (list) – concentrations of dose, drug concentration, and area under the curve

Returns:

dxdt (change of concentrations of dose, drug concentration, and area under the curve with t)

Return type:

list

calculate_PTH_concentration(t)

This function calculates the PTH concentration depending on time, basal PTH pulse and injected PTH pulse (if present).

Parameters:

t (float) – time variable

Returns:

PTH concentration

Return type:

float

calculate_PTH_concentration_from_basal(t, glandular_pulse)

This function calculates the PTH concentration based on the basal PTH pulse if no injected PTH is present. The PTH concentration depends on non-pulsatile (off-phase) or pulsatile (on-phase) share of the basal PTH pulse.

Parameters:

• t (float) – time variable

• glandular_pulse (int) – number of pulses for basal PTH

Returns:

PTH concentration

Return type:

float

calculate_PTH_concentration_from_basal_and_injected(t, glandular_pulse)

This function calculates the PTH concentration based on the basal and injected PTH pulse if injected PTH is present. The PTH concentration depends on non-pulsatile (off-phase) or pulsatile (on-phase) share of the basal PTH pulse and injected PTH pulse.

Parameters:

• t (float) – time variable

• glandular_pulse (int) – number of pulses for basal PTH

Returns:

PTH concentration

Return type:

float

calculate_activity_constants(cellular_activity, time)

This function calculates the activity constants of the cellular activity: basal activity (alpha_0), integrated activity (alpha_T), and cellular responsiveness (alpha_R). In the case of injected PTH, the activity constants are calculated for basal and injected PTH pulses. Otherwise, only the basal PTH pulse is considered.

Parameters:

• cellular_activity (list) – cellular activity (alpha(t) in original publication)

• time (list) – time variable

Returns:

basal activity, integrated activity, cellular responsiveness

Return type:

list

calculate_activity_constants_for_basal(basal_activity, integrated_activity_for_step_increase, cellular_activity, time)

This function calculates the activity constants integrated activity and cellular responsiveness for the basal PTH pulse if no injected PTH pulse is present. It first chooses one pulse of the cellular activity after adaptation for the calculation of the constants.

Parameters:

• basal_activity (float) – basal activity

• integrated_activity_for_step_increase (float) – integrated activity for a step increase of the same magnitude

• cellular_activity (list) – cellular activity depending on time

• time (list) – time variable

Returns:

integrated activity and cellular responsiveness

Return type:

list

calculate_activity_constants_for_basal_and_injection(basal_activity, integrated_activity_for_step_increase, cellular_activity, time)

This function calculates the activity constants integrated activity and cellular responsiveness for the basal and injected PTH pulse. It first chooses one pulse of the cellular activity after adaptation for both basal and injected for the calculation of the constants. The constants are calculated separately for basal and injected PTH pulses and then summed up.

Parameters:

• basal_activity (float) – basal activity

• integrated_activity_for_step_increase (float) – integrated activity for a step increase of the same magnitude

• cellular_activity (list) – cellular activity depending on time

• time (list) – time variable

Returns:

integrated activity and cellular responsiveness

Return type:

list

calculate_adaptation_time(stimulus_concentration)

This function calculates the adaptation time of the cellular activity to a stimulus depending on desensitised and resensitised receptors. It corresponds to the term tau_a in equation (11) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

adaptation time

Return type:

float

calculate_basal_activity()

This function calculates the basal activity of the cellular activity. It corresponds to the term alpha_0 in equation (9) in the Li & Goldbeter paper.

Returns:

basal activity

Return type:

float

calculate_cellular_activity()

This function calculates the cellular activity of the model based on receptor and complex concentrations. It corresponds to the term alpha(t) in equation (4) in the Li & Goldbeter paper.

Returns:

cellular activity and time

Return type:

list

calculate_cellular_responsiveness(integrated_activity, integrated_activity_for_step_increase, period)

This function calculates the cellular responsiveness of the cellular activity depending on the integrated activity. It corresponds to the term alpha_R in equation (31) in the Li & Goldbeter paper.

Parameters:

• integrated_activity (float) – integrated activity of the cellular activity

• integrated_activity_for_step_increase (float) – integrated activity for a step increase of the same magnitude

• period (float) – period of one PTH pulse

Returns:

cellular responsiveness

Return type:

float

calculate_contribution_of_receptor_desensitation(stimulus_concentration)

This function calculates the contribution of receptor desensitisation to the cellular activity. It corresponds to the term u in equation (12a) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

contribution of receptor desensitisation

Return type:

float

calculate_contribution_of_receptor_resensitisation(stimulus_concentration)

This function calculates the contribution of receptor resensitisation to the cellular activity. It corresponds to the term v in equation (12b) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

contribution of receptor resensitisation

Return type:

float

calculate_desensitised_receptors_after_adaptation(stimulus_concentration)

This function calculates the total fraction of desensitised receptors after adaptation to a stimulus. It depends on the contribution of receptor desensitisation and resensitisation. It corresponds to the term Ds in equation (8) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

desensitised receptors after adaptation

Return type:

float

calculate_difference_in_receptor_fraction(stimulus_concentration)

This function calculates the difference in receptor fraction for a stimulus concentration. It corresponds to the term Q in equation (20) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

difference in receptor fraction for minimum and maximum of PTH pulse

Return type:

float

calculate_difference_in_weights(stimulus_concentration)

This function calculates the difference in weights of active and desensitised receptors for a stimulus concentration. It corresponds to the term P in equation (10) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

difference in weights of active and desensitised receptors

Return type:

float

calculate_drug_PTH_pulse(load_case)

This function calculates the drug PTH pulse based on the load case. If there is a drug PTH injection, the ODE system (pharmacokinetics-pharmacodynamics model) is solved to determine the drug PTH pulse parameters for a square-wave pulse approximation (on/off-phase, pulse height).

Parameters:

load_case (Martonova_Load_Case) – Load case for the model

Returns:

None

Return type:

None

calculate_integrated_activity(chosen_activity_pulse, chosen_activity_pulse_time, basal_activity)

This function numerically calculates the integrated activity of one pulse of the cellular activity above baseline. It corresponds to the term alpha_T in equation (26) in the Li & Goldbeter paper.

Parameters:

• chosen_activity_pulse (list) – cellular activity of one pulse after adaptation

• chosen_activity_pulse_time (list) – time variable for the cellular activity

• basal_activity (float) – basal activity

Returns:

integrated activity

Return type:

float

calculate_integrated_activity_for_step_increase(stimulus_concentration)

This function calculates the integrated activity for a step increase of same magnitude in the cellular activity. It corresponds to the term alpha_Tstep in equation (28) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

integrated activity for step increase

Return type:

float

calculate_receptor_complex_concentrations()

This function calculates the concentrations of receptors and complexes over time by solving the ODE system for the receptor-ligand model.

Returns:

concentrations of active receptor, active complex, inactive complex, inactive receptor, and time

Return type:

list

calculate_weight_active_receptor(stimulus_concentration)

This function calculates the weight of active receptor in the cellular activity. It corresponds to the term a in equation (13a) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

weight of active receptor

Return type:

float

calculate_weight_desensitised_receptor(stimulus_concentration)

This function calculates the weight of desensitised receptor in the cellular activity. It corresponds to the term b in equation (13b) in the Li & Goldbeter paper.

Parameters:

stimulus_concentration (float) – stimulus/ ligand concentration

Returns:

weight of desensitised receptor

Return type:

float

receptor_ligand_model(t, x)

This function defines the receptor-ligand model based on an ODE system depending on the current PTH concentration. The system describes the change of concentrations of active receptor, active complex, and inactive complex over time. The concentration of inactive receptor is calculated based on the sum of the other concentrations. The receptor-ligand model corresponds to equation (2) in the Li & Goldbeter paper.

Parameters:

• t (float) – time variable

• x (list) – concentrations of active receptor, active complex, and inactive complex

Returns:

ODE system (change of concentrations of active receptor, active complex, and inactive complex)

Return type:

list

solve_for_activity()

This function calculates the time dependent cellular activity and the activity constants basal activity, cellular responsiveness, and integrated activity.

Returns:

cellular activity, time, basal activity, integrated activity, cellular responsiveness

Return type:

list

Modiz Model

class bone_models.bone_cell_population_models.models.modiz_model.Modiz_Model(load_case, model_type='cellular responsiveness', calibration_type='all')

Bases: Lemaire_Model

This class defines the bone cell population model semi-coupled with the 2-state receptor PTH model by Modiz et al., 2025. The model extends the bone cell population model by Lemaire et al., 2004, with activation of osteoblasts by PTH calculated by a 2-state receptor model with pulsatile PTH (Martonova et al., 2023). It is a subclass of the Lemaire model (see Lemaire_Model), inherits most methods, but modifies the calculate_PTH_activation_OB method. The Martonova model is used to calculate the activation by PTH separately for healthy and disease states.

Note

Source Publication: Modiz C., Castoldi N.M., Scheiner S., Martinez Reina J., Gallego J. L. C., Sansalone V., Martelli S., Pivonka P. (2025). Computational Simulations of Endocrine Bone Diseases Related to Pathological Glandular PTH Secretion Using a Multi-Scale Bone Cell Population Model Journal of Applied Mathematical Modelling (submitted)

The model is based on the following publications: Lemaire, V., Tobin, F. L., Greller, L. D., Cho, C. R., & Suva, L. J. (2004). Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. Journal of Theoretical Biology, 229(3), 293-309. DOI: 10.1016/j.jtbi.2004.03.023

Martonova D., Lavaill M., Forwood M.R., Robling A., Cooper D.M.L., Leyendecker S., Pivonka P. (2023). Effects of PTH glandular and external dosing patterns on bone cell activity using a two-state receptor model —Implications for bone disease progression and treatment PLOS ONE, 18(3), e0283544. DOI: 10.1371/journal.pone.0283544

Parameters:

• load_case (Modiz_Load_Cases) – load case for the model, include both load cases for Lemaire and Martonova model

• model_type (str) – type of the model (which activity constant is used) either ‘cellular responsiveness’ or ‘integrated activity’

• calibration_type (str) – type of calibration (alignment of activity constants and old activation using either all states or only healthy state), either ‘all’ or ‘only for healthy state’

• parameters (Modiz_Parameters) – instance of the Modiz_Parameters class

Raises:

• ValueError – If model_type is not ‘cellular responsiveness’ or ‘integrated activity’.

• ValueError – If calibration_type is not ‘all’ or ‘only for healthy state’.

Constructor method. Initialises parent class with respective load case and sets model type and calibration. Asserts that model type and calibration type are valid. Calculates the activity constants for healthy and disease states (in load case) using the Martonova model.

calculate_PTH_activation_OB(t)

Calculates the activation of osteoblasts by PTH. The calibration parameter is selected depending on the calibration type. The activation is calculated using either cellular responsiveness or integrated activity depending on the model type and the calibration parameter. Either healthy or disease activity constants are returned depending on the time and load case.

Parameters:

t (float) – time at which the activation is calculated, if None, the activation is calculated for the steady state

Returns:

activation of osteoblasts by PTH

Return type:

float

class bone_models.bone_cell_population_models.models.modiz_model.Reference_Lemaire_Model(load_case)

Bases: Lemaire_Model

This class is used to modify the Lemaire model to include a disease load case with multiplicative elevation. It is used a reference model for the calibration of the Modiz model.

Parameters:

load_case (Modiz_Load_Case) – load case for the model

Constructor method. Initialises parent class with respective load case.

calculate_PTH_concentration(t)

Calculates the PTH concentration. In a disease state (load case) the PTH concentration is elevated/ decreased by a respective multiplicative factor saved as load case parameter.

Parameters:

t (float) – time at which the PTH concentration is calculated

Returns:

PTH concentration

Return type:

float

bone_models.bone_cell_population_models.models.modiz_model.analyse_effect_of_different_pulse_characteristics(plot=True)

This function analyses the effect of different pulse characteristics on the integrated activity and cellular responsiveness. It calculates the integrated activity and cellular responsiveness for healthy and disease states using the Martonova model for different pulse characteristics. The pulse characteristics are varied by changing the on duration of the pulse on adn off phases. The integrated activity and cellular responsiveness are then plotted against the log of the ratio of the on and off duration of the pulse. The integrated activity and cellular responsiveness for the disease states are also plotted as a comparison.

Parameters:

plot (bool) – if True, the results are plotted

Returns:

Return type:

bone_models.bone_cell_population_models.models.modiz_model.calculate_elevation_parameter(disease)

Calculates the elevation/ reduction parameter for a disease state. The elevation parameter is the ratio of the minimum and maximum basal PTH pulse of the disease state to the minimum and maximum basal PTH pulse of the healthy state.

Parameters:

disease (Load_Case) – disease state

Returns:

elevation parameter

Return type:

float

bone_models.bone_cell_population_models.models.modiz_model.identify_calibration_parameters()

This function identifies the calibration (elevation/decrease) parameters for the Modiz model. It calculates integrated activity, cellular responsiveness and old activation for healthy and disease states using the Martonova model and the Lemaire model. The old activation of the Lemaire model is made comparable using the elevation/decrease parameter. It then performs an optimization to find the calibration parameters for cellular responsiveness and integrated activity using all states. The calibration parameters are printed. The calibration parameters are then saved as parameters for the Modiz model.

Print:

calibration parameters for cellular responsiveness and integrated activity

bone_models.bone_cell_population_models.models.modiz_model.identify_calibration_parameters_only_for_healthy_state()

This function identifies the calibration/alignment (elevation/decrease) parameters for the Modiz model. It calculates integrated activity, cellular responsiveness and old activation for only healthy state using the Martonova model and the Lemaire model. The parameter is then calculated for the healthy state only.

Print:

calibration parameters for cellular responsiveness and integrated activity

bone_models.bone_cell_population_models.models.modiz_model.objective_function(parameter, lemaire_activation_PTH, martonova_activation_PTH)

Objective function for the optimization to find the calibration parameters for cellular responsiveness and integrated activity.

Parameters:

• parameter (float) – calibration parameter

• lemaire_activation_PTH (list) – activation of osteoblasts by PTH calculated by the Lemaire model

• martonova_activation_PTH (list) – activation of osteoblasts by PTH calculated by the Martonova model

Returns:

error

Return type:

float

Pivonka Model

class bone_models.bone_cell_population_models.models.pivonka_model.Pivonka_Model(load_case)

Bases: Lemaire_Model

This class implements the bone cell population model by Lemaire et al. (2004) as a subclass of the Base_Model class.

Note

Source Publication: Pivonka P., Zimak J., Smith D.W., Gardiner B.S., Dunstan C.R., Sims N.A., Martin J.T., Mundy G.R. (2008). Model structure and control of bone remodeling: A theoretical study. Bone, 43(2), 249-263. DOI: 10.1016/j.bone.2008.03.025

Parameters:

• load_case (object) – load case for the model

• parameters (Pivonka_Parameters) – model parameters

• initial_guess_root (numpy.ndarray) – initial guess for the root-finding algorithm for steady-state

• steady_state (object) – steady state values of the model

• t (float) – time variable

• tspan (numpy.ndarray with start and end time) – time span for the ODE solver

• x (list) – state variables of the model

• solution (scipy.integrate._ivp.ivp.OdeResult) – solution of the ODE system

• initial_bone_volume_fraction (float) – initial bone volume fraction

• dOBpdt (float) – rate of change of precursor osteoblast cell concentration

• dOBadt (float) – rate of change of active osteoblast cell concentration

• dOCadt (float) – rate of change of osteoclast cell concentration

• dxdt (list) – rate of change of state variables

• bone_volume_fraction (list) – bone volume fraction over time

Constructor for the Lemaire_Model class.

Parameters:

load_case (object) – load case for the model

calculate_OPG_concentration(OBp, OBa, t)

Calculates the osteoprotegerin (OPG) concentration based on the osteoblasts’ production, PTH repression, maximum concentration, external injection and degradation rate. As described in the source publication, OPG can be produced by either precursor or active osteoblasts - determined by boolean variables.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

OPG concentration

Return type:

float

calculate_PTH_activation_OB(t)

Calculates the activation of osteoblasts by parathyroid hormone (PTH) based on the current PTH concentration.

Parameters:

t (float) – time variable

Returns:

activation of osteoblasts by PTH

Return type:

float

calculate_PTH_repression_OB(t)

Calculates the repression of osteoblasts by parathyroid hormone (PTH) based on the current PTH concentration.

Parameters:

t (float) – time variable

Returns:

repression of osteoblasts by PTH

Return type:

float

calculate_RANKL_RANK_concentration(OBp, OBa, t)

Calculates the RANKL-RANK complex concentration based on the RANKL concentration, RANK concentration and binding constant.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL-RANK complex concentration

Return type:

float

calculate_RANKL_activation_OCp(OBp, OBa, t)

Calculates the activation of precursor osteoclasts by RANKL based on the RANKL-RANK complex concentration and activation coefficient.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

activation of precursor osteoclasts by RANKL

Return type:

float

calculate_RANKL_concentration(OBp, OBa, t)

Calculates the RANKL concentration based on the effective carrying capacity, RANKL-RANK-OPG binding, degradation rate, intrinsic RANKL production and external injection of RANKL.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL concentration

Return type:

float

calculate_TGFb_activation_OBu(OCa, t)

Calculates the activation of uncommitted osteoblasts by TGF-beta based on the current TGF-beta concentration.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

activation of uncommitted osteoblasts by TGF-beta

Return type:

float

calculate_TGFb_activation_OCa(OCa, t)

Calculates the activation of active osteoclasts by TGF-beta based on the current TGF-beta concentration.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

activation of active osteoclasts by TGF-beta

Return type:

float

calculate_TGFb_concentration(OCa, t)

Calculates the TGF-beta concentration based on the osteoclastic resorption, external injection and degradation rate.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

TGF-beta concentration

Return type:

float

calculate_TGFb_repression_OBp(OCa, t)

Calculates the repression of precursor osteoblasts by TGF-beta based on the current TGF-beta concentration.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

repression of precursor osteoblasts by TGF-beta

Return type:

float

calculate_effective_carrying_capacity_RANKL(OBp, OBa, t)

Calculates the effective carrying capacity of RANKL based on the osteoblasts’ production, PTH activation and maximum RANKL per cell. As described in the source publication, RANKL can be produced by either precursor or active osteoblasts - determined by boolean variables.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

effective carrying capacity of RANKL

Return type:

float

calculate_external_injection_TGFb(t)

Calculates the external injection of TGF-beta (used and determined in load case scenarios).

Parameters:

t (float) – time variable

Returns:

concentration of external injection of TGF-beta

Return type:

float

Scheiner Model

class bone_models.bone_cell_population_models.models.scheiner_model.Scheiner_Model(load_case)

Bases: Pivonka_Model

This class implements the bone cell population model by Scheiner et al. (2013) as a subclass of the Pivonka_Model class. The main difference to the Pivonka model is the inclusion of mechanical effects on the precursor osteoblast cell concentration and RANKL concentration. The mutual influence between bone cells and mechanical effects is modeled as follows: Bone cells form and resorb bone matrix, which influences the strain energy density (mechanics model) felt by the bone cells (bone cell population model). The strain energy density in turn influences the proliferation rate of the osteoblast precursor cells and RANKL production. This again effects the cell concentrations and thus bone and vascular pore volume fraction.

Note

Source Publication: Scheiner, S., Pivonka, P., Hellmich, C. (2013). Coupling systems biology with multiscale mechanics, for computer simulations of bone remodeling. Computer Methods in Applied Mechanics and Engineering, 254, 181-196. DOI: 10.1016/j.cma.2012.10.015

Parameters:

• load_case (object) – load case for the model

• parameters (Scheiner_Parameters) – model parameters

• initial_guess_root (numpy.ndarray) – initial guess for the root-finding algorithm for steady-state

• steady_state (object) – steady state values of the model

• t (float) – time variable

• tspan (numpy.ndarray with start and end time) – time span for the ODE solver

• x (list) – state variables of the model

• solution (scipy.integrate._ivp.ivp.OdeResult) – solution of the ODE system

• initial_bone_volume_fraction (float) – initial bone volume fraction

• dOBpdt (float) – rate of change of precursor osteoblast cell concentration

• dOBadt (float) – rate of change of active osteoblast cell concentration

• dOCadt (float) – rate of change of osteoclast cell concentration

• dvascular_pore_fractiondt (float) – rate of change of vascular pore volume fraction

• dbone_volume_fractiondt (float) – rate of change of bone volume fraction

• dxdt (list) – rate of change of state variables

Constructor for the Scheiner_Model class as subclass of the Pivonka_Model class.

Parameters:

load_case (object) – load case for the model

apply_mechanical_effects(OBp, OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Applies the mechanical effects on the precursor osteoblast cell concentration. The mechanical effects depend on the strain energy density and proliferation rate of precursor osteoblasts. The latter is updated during the habitual loading phase.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• OCa (float) – active osteoclast cell concentration

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

rate of change of precursor osteoblast cell concentration due to mechanical effects

Return type:

float

bone_cell_population_model(x, t=None)

Calculates the system of ordinary differential equations for the bone cell population model, vascular pore volume fraction and bone volume fraction. This function is overwritten from the source model to add vascular pore volume fraction and bone volume fraction, that are necessary to solve in every time step.

Parameters:

• x (list) – state variables of the model

• t (float) – time variable

Returns:

rate of change of state variables

Return type:

list

calculate_RANKL_concentration(OBp, OBa, t)

Calculates the RANKL concentration based on the effective carrying capacity, RANKL-RANK-OPG binding, degradation rate, intrinsic RANKL production and external injection of RANKL. An additional RANKL production is added due to mechanical effects.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL concentration

Return type:

float

calculate_TGFb_concentration(OCa, t)

Calculates the TGF-beta concentration based on the osteoclastic resorption, external injection and degradation rate. Note: the resorption rate in this formula is included in the original model, but not coded.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

TGF-beta concentration

Return type:

float

calculate_hill_tensor_cylindrical_inclusion()

Calculates the fourth order Hill tensor for a cylindrical inclusion embedded in the bone matrix with a certain stiffness. The tensor is calculated by numerical integration using the stiffness tensor of the bone matrix (rewritten from matrix notation). The result is stored in the parameters object to avoid recalculation.

Returns:

Hill tensor for a cylindrical inclusion

Return type:

numpy.ndarray

calculate_macroscopic_stiffness_tensor(strain_concentration_tensor_bone_matrix, strain_concentration_tensor_vascular_pores, vascular_pore_fraction, bone_volume_fraction)

Calculates the macroscopic stiffness tensor depending on the strain concentration tensors, volume fractions and stiffness tensors for the both vascular pores and bone matrix. The calculation corresponds to Eq. (9) in the paper.

Parameters:

• strain_concentration_tensor_bone_matrix (numpy.ndarray) – strain concentration tensor for bone matrix

• strain_concentration_tensor_vascular_pores (numpy.ndarray) – strain concentration tensor for vascular pores

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

Returns:

macroscopic stiffness tensor

Return type:

numpy.ndarray

calculate_macroscopic_strain_tensor(strain_concentration_tensor_bone_matrix, strain_concentration_tensor_vascular_pores, vascular_pore_fraction, bone_volume_fraction, t)

Calculates the macroscopic strain tensor depending on the macroscopic stiffness tensor and the macroscopic stress vector. It is needed to calculate the microscopic strain tensor. The calculation corresponds to Eq. (13) in the paper.

Parameters:

• strain_concentration_tensor_bone_matrix (numpy.ndarray) – strain concentration tensor for bone matrix

• strain_concentration_tensor_vascular_pores (numpy.ndarray) – strain concentration tensor for vascular pores

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

macroscopic strain tensor

Return type:

numpy.ndarray

calculate_macroscopic_stress_vector(t)

Calculates the macroscopic stress vector depending on the load case scenario. The stress tensor is chosen depending on the time point (habitual loading phase or load case scenario). The tensor is the rewritten in vector form.

Parameters:

t (float) – time variable

Returns:

macroscopic stress vector

Return type:

numpy.ndarray

calculate_microscopic_strain_tensor(vascular_pore_fraction, bone_volume_fraction, t)

Calculates the microscopic strain tensor depending on the strain concentration tensors and the macroscopic strain tensor. The calculation corresponds to Eq. (14) in the paper.

Parameters:

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

microscopic strain tensor

Return type:

numpy.ndarray

calculate_strain_concentration_tensors(bone_volume_fraction)

Calculates the strain concentration tensors for the bone matrix and vascular pores depending on the hill tensor of the cylindrical inclusion, stiffness tensors for both phases and the current bone volume fraction. The calculation corresponds to Eq. (10) in the paper.

Parameters:

bone_volume_fraction (float) – bone volume fraction

Returns:

strain concentration tensor for bone matrix, strain concentration tensor for vascular pores

Return type:

numpy.ndarray, numpy.ndarray

calculate_strain_effect_on_OBp(OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Calculates the effect of strain on the proliferation rate of precursor osteoblasts. The effect depends on the time and strain energy density. For the steady-state, strain energy density is calculated once. During the load case scenario, the effect is updated based on the relation of current strain energy density to steady-state strain energy density. This relation also determines if RANKL production is set to 0 or computed depending on the strain energy density. The OBp proliferation rate is not updated during the load cse scenario.

Parameters:

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

effect of strain on the proliferation rate of precursor osteoblasts

Return type:

float

calculate_strain_energy_density(OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Calculates the microscopic strain energy density, experienced by the extravascular bone matrix, that drives the mechanoregulatory responses. It depends on the microscopic strain tensor (calculated in dependent functions) and the stiffness tensor of the bone matrix (fixed parameter). The calculation corresponds to Eq. (15) in the paper.

Parameters:

• vascular_pore_fraction (float) – vascular pore volume fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

microscopic strain energy density

Return type:

float

solve_bone_cell_population_model(tspan)

Solve the bone cell population model and volume fractions using the ODE system over a given time interval. The initial conditions are set to the steady-state values. This function is overwritten from the source model to add vascular pore volume fraction and bone volume fraction, that are necessary to solve in every time step.

Parameters:

tspan (numpy.ndarray with start and end time) – time span for the ODE solver

Returns:

solution of the ODE system

Return type:

scipy.integrate._ivp.ivp.OdeResult

stiffness_matrix_to_tensor()

Reshapes the stiffness matrix of the bone matrix into a 3x3x3x3 tensor.

Returns:

stiffness tensor of the bone matrix

Return type:

numpy.ndarray

Martinez-Reina Model

class bone_models.bone_cell_population_models.models.martinez_reina_model.Martinez_Reina_Model(load_case)

Bases: Scheiner_Model

This class implements the Martinez-Reina et al., 2019 model. It is a bone cell population model that includes precursor, active osteoblasts and active osteoclasts, the mechanical effects of bone remodelling in line with Scheiner et al., 2013, the mineralisation of bone with a queuing algorithm, and the effects of denosumab treatment and postmenopausal osteoporosis. The mechanical model is adjusted to account for trabecular bone and a Youngs modulus depending on the ash fraction. The bone volume resorbed and formed then agan depends on the bone cells stimulated by PMO, denosumab and mechnical stimulation.

Note

Source Publication: Martinez-Reina, J., & Pivonka, P. (2019). * Effects of long-term treatment of denosumab on bone mineral density: insights from an in-silico model of bone mineralization.* Bone, 125, 87-95. DOI: 10.1016/j.bone.2019.04.022

Parameters:

• load_case (Martinez_Reina_Load_Case) – load case of the model

• parameters (Martinez_Reina_Parameters) – model parameters

• initial_guess_root (numpy.ndarray) – initial guess for the root-finding algorithm for steady-state for OBp, OBa, OCa, vascular_pore_fraction, bone_volume_fraction

• steady_state (object) – steady state values of the model

• ageing_queue (numpy.ndarray) – queue for the mineralisation of bone

• denosumab_concentration_over_time (numpy.ndarray) – denosumab concentration over time

• time_for_denosumab (numpy.ndarray) – time for denosumab concentration

• bone_apparent_density (numpy.ndarray) – apparent density of bone over time

• bone_material_density (numpy.ndarray) – material density of bone over time

Constructor for the Scheiner_Model class as subclass of the Pivonka_Model class.

Parameters:

load_case (object) – load case for the model

calculate_RANKL_concentration(OBp, OBa, t)

Calculates the RANKL concentration based on the effective carrying capacity, RANKL-RANK-OPG binding, degradation rate, intrinsic RANKL production and external injection of RANKL (PMO) and denosumab effect. Postmenopausal osteoporosis is simulated via an external injection of RANKL (see Eq. (18) in the source paper). Denosumab effect is included based on the current concentration in the compartment, binding constant of denosumab and RANKL and an accessibility factor.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL concentration

Return type:

float

calculate_ash_fraction(OBa, OCa, bone_volume_fraction, t)

Calculate the ash fraction based on the average mineral content, the volume fraction of mineral and organic material and the respective densities. The volume fraction of mineral is updated before the calculation.

: param OBa: active osteoblast cell concentration : type OBa: float : param OCa: active osteoclast cell concentration : type OCa: float : param bone_volume_fraction: bone volume fraction : type bone_volume_fraction: float : param t: time variable : type t: float

: return: ash fraction : rtype: float

calculate_average_mineral_content(OBa, OCa, bone_volume_fraction, t)

Calculate the average mineral content based on the mineralisation law and the mineralisation queue. Each element of the queue is multiplied by the mineral content determined by mineralisation law and summed up. The last element is assigned the maximum mineral content. The average mineral content is then calculated by dividing the sum by the bone volume fraction.

: param OBa: active osteoblast cell concentration : type OBa: float : param OCa: active osteoclast cell concentration : type OCa: float : param bone_volume_fraction: bone volume fraction : type bone_volume_fraction: float : param t: time variable : type t: float

: return: average mineral content : rtype: float

calculate_compliance_matrix(OBa, OCa, bone_volume_fraction, t)

Calculate the compliance matrix based on the current Young’s modulus and Poisson’s ratio. The Youngs modulus is updated before the calculation.

Parameters:

• OBa (float) – active osteoblast cell concentration

• OCa (float) – active osteoclast cell concentration

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

compliance matrix

Return type:

numpy.ndarray

calculate_denosumab_concentration(t)

Calculates the denosumab concentration based on the current time and the pharmacokinetics-pharmacodynamics model. If the time is outside of the treatment period, the concentration is zero. Otherwise, the concentration is calculated using the solution of the PK-PD model after translating the time back to the solution time interval. The PK-PD model returns the concentration in ng/ml, which is then converted to pmol/L using the molar mass of denosumab.

Parameters:

t (float) – time variable

Returns:

denosumab concentration

Return type:

float

calculate_external_injection_RANKL(t)

Calculate the external injection of RANKL based on the time variable for simulation of PMO. The injection is zero if the time is outside of the PMO simulation period. Otherwise, the injection is calculated based on the increase in RANKL, the reduction factor, the characteristic time and the time variable. The RANKL increase due to PMO is assumed to start at a fixed value and decrease over time.

Parameters:

t (float) – time variable

Returns:

external injection of RANKL (increase due to PMO)

Return type:

float

calculate_mineralisation_law(t)

Calculate the mineralisation law based on the time variable. In the lag time, the mineral content is zero. In the primary phase, the mineral content increases linearly from zero to the primary mineral content. In the secondary phase, the mineral content increases exponentially to the maximum mineral content. The mineral content is based on the age of the patch in the queue.

Parameters:

t (float) – time variable

Returns:

mineral content

Return type:

float

calculate_strain_energy_density(OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Calculate the strain energy density based on macroscopic stress vector, compliance matrix, strain matrix and stiffness tensor.

Parameters:

• OBa (float) – active osteoblast cell concentration

• OCa (float) – active osteoclast cell concentration

• vascular_pore_fraction (float) – vascular pore fraction

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

strain energy density

Return type:

float

calculate_youngs_modulus(OBa, OCa, bone_volume_fraction, t)

Calculate the Young’s modulus based on the bone volume fraction and ash fraction. The ash fraction is updated before the calculation.

Parameters:

• OBa (float) – active osteoblast cell concentration

• OCa (float) – active osteoclast cell concentration

• bone_volume_fraction (float) – bone volume fraction

• t (float) – time variable

Returns:

Young’s modulus

Return type:

float

initialise_ageing_queue(OBa, OCa, bone_volume_fraction)

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

pharmacokinetics_pharmocodynamics_denosumab(t, concentration)

The PK-PD model for denosumab concentration over time. The ODE is based on the absorption rate, elimination rate, volume of the central compartment, bioavailability, Michaelis-Menten constant, maximum volume and injected dose.

Parameters:

• t (float) – time variable

• concentration (float) – injected denosumab concentration

Returns:

derivative of the denosumab concentration

Return type:

float

solve_bone_cell_population_model(tspan)

Solve the bone cell population model and volume fractions using the ODE system over a given time interval. The initial conditions are set to the steady-state values. The queuing algorithm should only be updated every day, which prohibits the evaluation in each time step of the solver. Thus, the bone cell population model is solved for each interval [day, day+1] and the queuing algorithm is updated at the end of each interval with the current formed and resorbed bone volume fraction. Additionally, the bone material/ apparent density are calculated. The return array includes the time vector, the concentrations of OBp, OBa, OCa, the vascular pore fraction and the bone volume fraction.

Parameters:

tspan (numpy.ndarray with start and end time) – time span for the ODE solver

Returns:

solution of the ODE system [time, OBp, OBa, OCa, vascular_pore_fraction, bone_volume_fraction]

Return type:

list of arrays

solve_for_denosumab_concentration()

Solve the PK-PD model for denosumab concentration over time. The initial concentration is set to zero and the ODe is solved over the treatment period.

Returns:

None

update_ageing_queue(OBa, OCa, bone_volume_fraction)

Update the ageing queue based on the current resorbed and formed bone volume fraction. The queue is updated by shifting all elements by one position, resorb volume based on the elements size and adding the new formed bone volume fraction at the beginning. The last element of the queue stores the volume needed for all the elements to sum up to the bone volume fraction.

Parameters:

• OBa (float) – active osteoblast cell concentration

• OCa (float) – active osteoclast cell concentration

• bone_volume_fraction (float) – bone volume fraction

Returns:

None

Lerebours Model

class bone_models.bone_cell_population_models.models.lerebours_model.Lerebours_Model(load_case, porosity, specific_surface_multiplier=1)

Bases: Scheiner_Model

Implements the Lerebours mechanobiological model for bone cell population dynamics.

This model extends the Scheiner framework by incorporating porosity-dependent feedback loops. It couples bone cell populations (osteoblasts/osteoclasts) to the bone volume fraction via the specific surface area, allowing for simulations of site-specific bone loss and remodeling.

Note

Source Publication: Lerebours, C., Buenzli, P. R., Scheiner, S., & Pivonka, P. (2016). A multiscale mechanobiological model of bone remodeling predicts site-specific bone loss in the femur during osteoporosis and mechanical disuse. Biomechanics and Modeling in Mechanobiology, 15(1), 43-67. DOI: 10.1007/s10237-015-0705-x

Initializes the model with local geometry and loading conditions.

Parameters:

• load_case (Load_Case) – Object containing mechanical loading parameters.

• porosity (float) – Initial bone porosity (0 to 1), used to derive the specific surface area for cell recruitment.

• specific_surface_multiplier (float, optional) – Adjusts the surface availability for cell attachment. Defaults to 1.

apply_mechanical_effects(OBp, OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Computes the mechanical stimulus-driven change in osteoblast precursor (OBp) proliferation.

This method applies a piecewise mechanotransduction function to the baseline proliferation rate. The effect is determined by the strain energy density and the biomechanical transduction strength defined in the model parameters.

Parameters:

• OBp (float) – Precursor osteoblast concentration.

• OBa (float) – Active osteoblast concentration.

• OCa (float) – Active osteoclast concentration.

• vascular_pore_fraction (float) – Vascular pore volume fraction.

• bone_volume_fraction (float) – Bone volume fraction.

• t (float) – Time variable.

Returns:

Rate of change of OBp concentration due to mechanical effects.

Return type:

float

bone_cell_population_model(x, t=None)

Defines the ODE system for bone cell populations and volume fractions.

This function supports two modes: 1. Steady State (t=None): Calculates residuals for uncommitted and precursor cells given fixed active cell densities. 2. Transient (t=float): Calculates rates of change for precursors, active cells, and volume fractions.

Parameters:

• x (list) – State vector. If t is None: [OBu, OBp, OCu, OCp, vascular_pore_fraction, bone_volume_fraction]. If t is float: [OBp, OBa, OCp, OCa, vascular_pore_fraction, bone_volume_fraction].

• t (float or None) – Time variable; if None, the model assumes a steady-state calculation.

Returns:

List of rates of change [dOBp/dt, dOBa/dt, dOCp/dt, dOCa/dt, dVPF/dt, dBVF/dt].

Return type:

list

calculate_MCSF_activation_OCu()

Calculates the MCSF activation for uncommitted osteoclast cells (OCu) based on the MCSF concentration (constant parameter) and activation coefficient (Hill function).

Returns:

MCSF activation for uncommitted osteoclast cells (OCu)

Return type:

float

calculate_OPG_concentration(OBp, OBa, t)

Calculates the OPG concentration based on the effective carrying capacity of OPG, OPG production, calibration parameter, PTH repression and external injection of OPG.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

OPG concentration

Return type:

float

calculate_PTH_concentration(t)

Calculates the PTH concentration based on intrinsic PTH production and external injection of PTH.

Parameters:

t (float) – time variable

Returns:

PTH concentration

Return type:

float

calculate_RANKL_activation_OCp(OBp, OBa, t)

Calculates the RANKL activation for osteoclast precursor cells (OCp) based on the RANKL concentration. It calls the function calculate_RANKL_concentration.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL activation for osteoclast precursor cells (OCp)

Return type:

float

calculate_RANKL_activation_OCu(OBp, OBa, t)

Calculates the RANKL-mediated activation of uncommitted osteoclasts (OCu). This implementation assumes the activation rate for uncommitted cells is identical to that of osteoclast precursors (OCp).

Parameters:

• OBp (float) – Concentration of precursor osteoblasts.

• OBa (float) – Concentration of active osteoblasts.

• t (float) – Time variable.

Returns:

RANKL activation factor for OCu.

Return type:

float

calculate_RANKL_concentration(OBp, OBa, t)

Calculates the RANKL concentration based on the effective carrying capacity, RANKL-RANK-OPG binding, degradation rate, intrinsic RANKL production and external injection of RANKL. An additional RANKL production is added due to mechanical effects.

Parameters:

• OBp (float) – precursor osteoblast cell concentration

• OBa (float) – active osteoblast cell concentration

• t (float) – time variable

Returns:

RANKL concentration

Return type:

float

calculate_TGFb_concentration(OCa, t)

Calculates the TGFb concentration based on the stored TGFb content in bone volume, active osteoclasts, resorption rate, degradation rate and calibration parameter.

Parameters:

• OCa (float) – active osteoclast cell concentration

• t (float) – time variable

Returns:

TGFb concentration

Return type:

float

calculate_steady_state(porosity)

Calculates the steady state for the bone cell population model based on the porosity. It determines the active osteoclasts and osteoblasts based on the turnover, and then solves the bone cell population model for steady state of uncommitted and precursor cells. The steady state concentrations are stored as parameters (uncommitted stays constant for all further calculation).

Parameters:

porosity (float) – Porosity of the bone

Returns:

None

calculate_strain_effect_on_OBp(OBa, OCa, vascular_pore_fraction, bone_volume_fraction, t)

Calculates the normalized mechanical stimulus based on strain energy density (SED). In steady-state (t=None), this method initializes the baseline SED. During loading scenarios, it computes the relative deviation from this baseline to determine the mechanical effect on OBp proliferation and adapt RANKL production rates.

Parameters:

• OBa (float) – Active osteoblast concentration.

• OCa (float) – Active osteoclast concentration.

• vascular_pore_fraction (float) – Vascular pore volume fraction.

• bone_volume_fraction (float) – Bone volume fraction.

• t (float or None) – Time variable; if None or before start_time, returns 0 (steady-state).

Returns:

Normalized strain effect (relative change in SED).

Return type:

float

calculate_turnover(porosity)

Calculates the turnover based on porosity, calibration factor and specific surface. The calibration factor was identified by fitting the SV curve to turnover datapoints.

Parameters:

porosity (float) – Porosity of the bone

Returns:

Turnover rate

Return type:

float

set_macroscopic_stress_tensor(stress_xx, stress_xy, stress_xz, steady_state=False)

Sets the macroscopic stress tensor and updates the model or load case parameters.

Note

Coordinate System Convention: To maintain compatibility between the Scheiner and Lerebours frameworks, the input stress_xx is mapped to the internal z-direction (tensor index [2,2]). This discrepancy is essential for the correct operation of the spatial model component.

Parameters:

• stress_xx (float) – Axial stress (mapped to the internal z-direction).

• stress_xy (float) – Shear stress component in the xy plane.

• stress_xz (float) – Shear stress component in the xz plane.

• steady_state (bool) – If True, updates the normal loading parameters; otherwise, updates the current load case.

Returns:

None

solve_bone_cell_population_model(tspan, porosity, initial_conditions=None)

Solve the bone cell population model and volume fractions using the ODE system over a given time interval. The initial conditions are set to the steady-state values which are calculated if initial conditions are None. This function is overwritten from the source model to add vascular pore volume fraction and bone volume fraction, that are necessary to solve in every time step.

Parameters:

• tspan (numpy.ndarray with start and end time) – time span for the ODE solver

• porosity (float) – Porosity of the bone ([0,1]), used to calculate specific surface and steady state

• initial_conditions (list or None) – Initial conditions for the ODE system, if None, steady state is calculated

Returns:

solution of the ODE system

Return type:

scipy.integrate._ivp.ivp.OdeResult

specific_surface(porosity)

This function calculates the specific surface of bone based on the porosity. The specific surface multiplier is used for the boundaries in Modiz et al. (irrelevant for this model), default is 1.

Parameters:

porosity (float) – Porosity of the bone

Returns:

Specific surface of bone

Return type:

float