Patent Application
20210357548
The present teachings are generally related to methods and systems for performing in-silico modeling of bone remodeling, and in particular, methods and systems for in-silico modeling of treatment and progression of osteoporosis.
Bone remodeling is a process by which bone strength and homeostatic function of inorganic components are maintained. The bone remodeling involves replacing portions of the old bone with newly formed bone matrix. Osteoclasts are responsible for bone resorption and osteoblasts are responsible for secreting new bone matrix in a process called ossification. Osteocytes are essentially osteoblasts surrounded by the bone matrix.
Osteoporosis is a disease in which bone density (quantitated as bone mineral density) deceases and results in a porous bone matrix that is more susceptible to fracture. Osteoporosis can result from an increase in bone matrix resorption and/or a decrease in the formation of new bone matrix. A variety of factors can affect the progression of osteoporosis. Some examples of such factors include the age of onset of menopause, age, medication, physical activity, among others. For example, bone resorption rates are much higher in post-menopausal women due to estrogen deficiency caused by menopause. Common treatments of osteoporosis include the administration of drugs that increase bone mineral density. Some examples of such drugs include bisphosphonates, RANKL inhibitors, and selective estrogen receptor modulators (SERMs). Osteopenia is characterized by reduced bone mass of lesser severity than osteoporosis, and can be precursor of osteoporosis.
It is important to know the effect of certain physiological conditions as well as certain therapeutic protocols on the treatment and progression of osteoporosis. Clinical trials can be performed to obtain such information. But conducting a clinical trial is usually expensive and time-consuming.
Accordingly, there is a need for cost-effective, predictive models that can be used to obtain information regarding bone remodeling, and in particular, the effects of various physiological conditions and therapeutic protocols on the treatment and progression of osteoporosis.
In one aspect, a method for optimizing a therapy for osteoporosis is disclosed, which comprises (a) generating in silico a plurality of virtual patients, wherein each of the virtual patients comprises a mathematical construct for modeling progression of osteoporosis via simulation of bone remodeling, where the mathematical construct comprises a plurality of dynamic mathematical relations defining time-dependent evolution of at least one of a cell density variable or concentration variable (e.g., cell population) associated with any of pre-osteoblasts, osteoblasts, preosteoclasts, osteoclasts, osteocytes, and a bone resorption signal, sclerostin, estrogen, bone density and bone mineral content and employs a processor to determine time-variation of at least one of bone mineral density and bone mineral fraction based on at least a portion of said time-dependent variables, (b) applying a simulated therapy to said plurality of virtual patients, and (c) using said virtual patients to determine an effect of said simulated therapy on progression of osteoporosis.
In some embodiments, the method can further include: (d) adjusting the simulated therapy by adjusting at least one of the parameters of at least one of the dynamic mathematical relations, (e) applying said adjusted simulated therapy to said virtual patients, and (f) determining an effect of said adjusted simulated therapy on the progression of osteoporosis.
In some embodiments, the above steps d)-f) can be repeated so as to optimize the simulated therapy.
In some embodiments, at least one parameter of at least one of the dynamic mathematical relations is determined based on physiological data collected from a population of actual patients. By way of example, the physiological data can include any of measured bone mineral density, and measured serum concentration of one or more bone turnover markers. In some embodiments, one or more bone turnover markers comprise any of N-terminal propeptide of type I procollagen (PINP), bone-specific alkaline phosphatase (BSAP), and C-terminal telopeptide (CTX).
In some embodiments, the simulated therapy can include the virtual administration of at least one medication to the virtual patients. In some such embodiments, the simulated therapy includes the virtual administration of a plurality of medications (e.g., two or more medications) to the virtual patients. In some such embodiments, the pharmacokinetic data associated with the one or more medications is employed to simulate the effect of those medications (e.g., individually or in combination) on the progression of osteoporosis. Further, in some such embodiments, the effect of the simulated therapy on the progression of osteoporosis can be determined based on a particular dosage and/or regimen employed for the administration of one or more medications to the virtual patients.
In some embodiments, the simulated therapy includes adjusting a dosage of at least one medication (or multiple medications) that is virtually administered to the virtual patients. For example, such adjustment of the medication can include adjusting any of the dose and/or the daily schedule for the virtual administration of the medication(s) to the virtual patients.
Further, in some embodiments, at least one of the virtual patients is configured to simulate an effect of a disease on the progression of the osteoporosis. By way of example, the disease can be multiple myeloma. In some such embodiments, at least one of the virtual patients simulates the effect of multiple myeloma on the progression of osteoporosis by adjusting at least one of the following rates: a birth rate, a differentiation rate or cell death rate of at least one of pre-osteoblasts, pre-osteoclasts, osteoblasts, osteoclasts, and osteocytes. In some such embodiments, the step of adjusting at least one of the rates includes defining a mathematical relation indicative of time variation of the rate after onset of the disease.
In another embodiment, the effect of glucocorticoid-induced osteoporosis on the progression of osteoporosis can be simulated by adjusting at least one of the following rates: a birth rate, a differentiation rate or cell death rate of at least one of pre-osteoblasts, pre-osteoclasts, osteoblasts, osteoclasts, and osteocytes.
In some embodiments, at least one of the variables can be adjusted to simulate an effect of a change in a physiological condition on the progression of osteoporosis. The physiological condition can be, for example, the onset of menopause. In some embodiments, the estrogen concentration can be used as the interface variable for simulating the effect of menopause on the progression of osteoporosis. In some such embodiments, the simulation of the effect of menopause on the progression of osteoporosis can be performed by using a time-dependent relation defining post-menopausal decline of estrogen as a function of time.
In some embodiments, the physiological condition is caused by aging. In some embodiments, the variable associated with the age-dependent physiological condition comprises the serum concentration of sclerostin. By way of example, the simulation of the effect of aging can be achieved by providing a relation that defines a time-dependent increase in the serum concentration of the sclerostin after a selected age threshold.
In some embodiments, the method can further include adjusting at least one of the dynamic variables to simulate the effect of application of a therapy on the progression of osteoporosis. In one embodiment, the therapy includes the virtual administration of a RANKL inhibitor to an individual and the step of adjusting at least one variable includes defining a mathematical relation indicative of time variation of any of the pre-osteoclast cell density, the osteoclast cell density and the bone mineral content as a function of the virtual administration of the RANKL inhibitor.
Further, in some embodiments, the therapy includes the virtual administration of a sclerostin inhibitor and the step of adjusting at least one of the dynamic variables includes defining a mathematical relation that is indicative of time variation of the sclerostin concentration as a function of the virtual administration of the Sclerostin inhibitor.
In some embodiments, the therapy includes the virtual administration of bisphosphonates and the step of adjusting at least one variable comprises defining a mathematical relation indicative of time variation of the osteoclasts cell density as a function of enhanced apoptosis rate due to the virtual administration of bisphosphonates. By way of example, the bisphosphonates can be, without limitation, any of alendronate, risedronate, ibandronate, and zoledronate.
In some embodiments, the method includes adjusting at least one of the dynamic variables to simulate the effect of a physical activity and/or microgravity on any of treatment and progression of osteoporosis. In some such embodiments, the at least one variable associated with the physical activity includes a concentration of the sclerostin and an in-silico model according to the present teachings provides a mathematical relation indicating a temporal change in the concentration of sclerostin as a function of time due to physical activity, as discussed in more detail below. In a related aspect, a method for optimizing a therapy for osteoporosis is disclosed, which includes: (a) generating in silico at least one virtual patient as a mathematical construct for modeling age-related progression of osteoporosis via simulation of bone remodeling, said mathematical construct comprising a plurality of dynamic mathematical relations defining time-dependent evolution of at least one of a cell density variable or concentration variable associated with any of pre-osteoblasts, osteoblasts, preosteoclasts, osteoclasts, osteocytes, a bone resorption signal, sclerostin, estrogen, bone density and bone mineral content and employs a processor to determine time-variation of bone mineral density based on at least a portion of said time-dependent variables, (b) applying a simulated therapy to said at least one virtual patient, and (c) using said at least one virtual patient to determine an effect of said simulated therapy on progression of osteoporosis. In some embodiments, the at least one virtual patient includes a plurality of virtual patients.
In some embodiments, the mathematical construct employs a relation indicative of variation of any of estrogen and sclerostin as a function of age to determine the effect of aging on at least one of said variables. By way of example, such a relation can be as follows:
where ae denotes the age at the onset of estrogen decline and τe is a characteristic time scale.
In some embodiments, the above method can further include adjusting the simulated therapy by adjusting at least one of said parameters of at least one of the dynamic mathematical relations, applying the adjusted simulated therapy to the virtual patients, and determining an effect of the adjusted simulated therapy on the progression of osteoporosis. This process can be iterated so as to obtain an optimal therapy for administration to the patient.
In some embodiments, at least one parameter of at least one of the dynamic mathematical relations is determined based on physiological data collected from a population of actual patients.
In some embodiments, the physiological data include any of measured bone mineral density and measured serum concentration of one or more bone turnover markers. Some examples of such markers include, without limitation, any of N-terminal propeptide of type I procollagen (PINP), bone-specific alkaline phosphatase (B SAP), and C-terminal telopeptide (CTX).
In some embodiments, the simulated therapy includes the virtual administration of at least one medication to the virtual patients. In some such embodiments, the simulated therapy includes the virtual administration of a combination of two or more medications to the virtual patients. In some embodiments, pharmacokinetic data associated with at least one medication is employed for determining the effect of a simulated therapy on the progression of osteoporosis.
In some embodiments, the adjustment of the simulated therapy includes adjusting a dose of at least one medication administered to the virtual patients. For example, the dose and/or the regimen of the administration of the medication to the virtual patients can be modified and the effect of the modification on the progression of osteoporosis can be simulated.
In a related aspect, a computer system for modeling osteoporosis is disclosed, which comprises a simulation module using at least one computer processor for generating in silico a mathematical construct comprising a plurality of dynamic mathematical relations for simulating bone remodeling, where the dynamic mathematical relations define time-dependent evolution of at least one of a density variable or a concentration variable associated with any of pre-osteoblasts, osteoblasts, pre-osteoclasts, osteoclasts, osteocytes, a bone resorption signal, sclerostin, estrogen, bone density and bone mineral content, said simulation module further being configured to determine time-variation of bone mineral density based on at least a portion of the time-dependent variables. The system can further include an avatar module for storing patient-specific parameters and being in communication with the simulation module to provide the patient-specific parameters thereto so as to enable the simulation module to simulate the effect of any of aging, physical activity, a physiological condition, a therapeutic treatment, a disease, a dietary or substance intake on progression of osteoporosis, and a fitting module configured to derive said patient-specific parameters by fitting one or more observables to the clinical data. By way of example, the clinical data can be based on data collected from a population of patients.
In some embodiments, the observables can include any of bone mineral density and bone turnover markers.
In some embodiments, the simulation module is configured to simulate the effect of multiple myeloma on the progression of osteoporosis. For example, the simulation module can simulate the effect of multiple myeloma on the progression of osteoporosis by adjusting a reference differentiation rate of the pre-osteoclasts. The adjustment of the reference differentiation rate includes defining a mathematical relation indicative of time variation of the reference differentiation rate after the onset of the disease.
In some embodiments, the simulation module is configured to simulate the effect of glucocorticoid-induced osteoporosis on the progression of osteoporosis by adjusting a reference differentiation rate of the pre-osteoclasts. In some embodiments, the physiological condition is onset of menopause, and the associated dynamic variable for simulating the effect of the onset of menopause on the progression of osteoporosis is the estrogen concentration. In some such embodiments, the simulation module employs a time-dependent relation defining post-menopausal decline of estrogen as a function of time to simulate the effect of menopause on the progression of osteoporosis.
In some embodiments, the simulation module is configured to adjust a concentration of sclerostin to simulate the effect of aging on the progression of osteoporosis. In some such embodiments, the simulation module employs a relation defining time-dependent increase in the serum concentration of the sclerostin after a selected age threshold to simulate the effect of aging on the progression of osteoporosis.
In some embodiments, the simulation module is configured to simulate the effect of the administration of a RANKL inhibitor on the progression of osteoporosis. In some such embodiments, the simulation module performs the simulation by defining a mathematical relation that is indicative of time variation of any of the pre-osteoclast cell density, the osteoclast cell density and the bone mineral content as a function of the administration of the RANKL inhibitor.
In some embodiments, the therapeutic treatment comprises the administration of a sclerostin inhibitor and the step of adjusting at least one variable comprises defining a mathematical relation indicative of time variation of the sclerostin concentration as a function of the virtual administration of the sclerostin inhibitor.
In some embodiments, the therapeutic treatment includes the administration of bisphosphonates and the step of adjusting at least one variable includes defining a mathematical relation indicative of time variation of the osteoclasts cell density as a function of enhanced apoptosis rate due to the virtual administration of bisphosphonates.
In a related aspect, a computer system for modeling osteoporosis is disclosed, which comprises a fitting module, an avatar module, and a simulation module. The fitting module is configured to receive any of patient-level and population-level clinical data. The fitting module also includes at least one cost function for determining whether a set of calibration parameters associated with the avatar module are acceptable, as discussed in more detail below. The fitting module is further configured to define fit parameters and parameter bounds and communicate this information to the avatar module. The avatar module is in turn configured to pass this information to the simulation module, which is configured to employ the dynamic mathematical relations of an in-silico model according to the present teachings to compute time evolution of one or more variables associated with bone remodeling and pass this information back to the avatar module.
The avatar module employs this information to compute one or more avatar-specific quantities, such as time-dependence of the bone turnover markers. The avatar module compiles the information received from the simulation module and the information it itself computes and transmits the compiled information to the fitting module.
The fitting module employs the compiled information to compute a value of the cost function and compares a reduction in the value of the cost function in each iteration with a predefined threshold. If the value is less than the threshold, the avatar calibration is complete and the fitting module can invoke the calibrated avatar module for obtaining time evolution of parameters of interest. If the value is greater than the threshold, the fitting modules adjusts the avatar parameters and passes the new parameters to the avatar module. This process can be reiterated until the value of the cost function is within an acceptable range.
The present teachings are generally directed to methods and systems for in-silico modeling of bone remodeling, and in particular the progression and/or treatment of osteoporosis, which is a medical condition characterized by a decrease in bone mass with decreased density and enlargement of bone spaces producing porosity and fragility.
Various terms are used herein in accordance with their ordinary meanings in the art. The term “about” as used herein denotes a variation of at most 10% around a numerical value and the term “substantially” as used herein denotes a deviation of at most 10% with respect to a complete state and/or condition. The term “osteoporosis” as used herein is intended to encompass any bone loss disease, and illness condition including, without limitation, osteopenia and more advanced stages of bone loss that are generally referred to in medicine as osteoporosis.
As discussed in more detail below, in some embodiments, the modeling of osteoporosis employs a mathematical construct implemented on a computing device, where in some embodiments, the mathematical construct includes a plurality of dynamic mathematical relations that define time-dependent evolution of at least one of a cell density variable or a cell population variable associated with any of pre-osteoblasts, osteoblasts, osteocytes, or a variable indicative of a bone resorption signal, sclerostin and/or estrogen concentration, bone density and bone mineral content.
In some embodiments, the mathematical construct includes only those components of the bone mineral metabolism that is relevant to disease progression and intervention. In particular, in many embodiments, the present teachings employ key modulators that are especially important in the context of osteoporosis (e.g., sex hormones and medication).
In many embodiments, concentrations and densities of various components of the model (e.g., osteoblasts, osteoclasts, etc.) are normalized to steady state healthy adult values in order to provide a non-dimensional model. In some embodiments, a number of time constants (e.g., equilibrium rates encoded in the degradation/apoptosis rates) can be inferred from known values in the literature. Subsequently, the gain rates can be determined self-consistently by imposing a unit steady state that represents the average healthy adult condition.
For self-consistency, in some embodiments, all components that do have a regulatory effect working in both directions are assigned threshold concentrations that are in the range of the steady-state value within an order of magnitude.
Further, in many embodiments, the effects of external factors on the RANK-RANKL-OPG pathway are taken into account by their effective action on the respective cell populations.
As discussed in more detail below, the following features, among others, are incorporated in an in-silica model according to the present teachings for bone mineralization/osteoporosis.
RANK-RANKL-OPG pathway is central to balancing bone formation and bone resorption. RANKL, is a ligand that can bind to RANK receptor of osteoclasts to stimulate them, thereby causing bone resorption. OPG (osteoprotegerin) is secreted by osteoblasts and binds to RANKL to prevent it from binding to RANK receptor, thereby inhibiting excessive bone resorption.
In bone matrix, several growth factors are present in a latent form. Examples of such growth factors include TGFβ, BMF2, PDGF and IGFs. These growth factors are deposited by the osteoblasts upon bone formation and are released by osteoclasts upon bone resorption. It is likely that they act on osteoblast progenitor populations. Estrogen inhibits bone resorption by directly inducing apoptosis of the bone-resorbing osteoclasts, among other processes. For example, estrogen is also known to regulate sclerostin, pre-osteoclast differentiation, and osteoblast apoptosis.
Sclerostin is a small protein expressed by the SOST gene in osteocytes. When sclerostin binds to its receptors on the cell surface of osteoblasts, a downstream cascade of intracellular signaling is initiated to inhibit bone formation. Circulating levels of sclerostin are regulated by estrogen. From the age 20 to 80, sclerostin serum levels can increase up to about fourfold. However, sclerostin mRNA levels in bone show no significant age dependence. In postmenopausal women, the serum level of sclerostin increases with a concomitant decrease in bone density. In postmenopausal women, the order of magnitude of serum sclerostin levels is about 0.3 ng/ml.
As discussed in more detail below, the model (and associated in silico system) can provide modules for simulating the effect of a disease and/or a therapeutic treatment on the bone mineralization. By way of example,
In the embodiments discussed below, the following conventions and abbreviations are employed.
Superscripts indicate components; e.g., ρB denotes osteoblast density. Regulatory threshold concentrations are denoted by the same symbol as the dynamic variable with a subscript indicating on which component the regulation acts; e.g., sB is the threshold sclerostin concentration for the regulation of osteoblasts (B).
In the discussion below, the following units are employed:
In some embodiments, systems and methods according to the present teachings for modeling bone mineralization employ a plurality of dynamic mathematical relations for determining the time evolution of any of concentration and density of various components of the model.
In some embodiments, bone remodeling can be modeled in silica by defining a plurality of dynamic mathematical relations that incorporate the following information regarding the transitions of pre-osteoblasts to osteoblasts and osteoblasts to osteocytes: (a) pre-osteoblast to osteoblast differentiation is upregulated by the resorption signal and downregulated by sclerostin; (b) osteoblast apoptosis is downregulated by estrogen; (c) osteoblast to osteocyte transition is downregulated by sclerostin; (d) osteocytes are exclusively generated via transition from osteoblasts and undergo apoptosis at a constant rate, and (e) osteoblast to osteocyte transition is downregulated by sclerostin. More specifically, in an embodiment of an in-silico model according to the present teachings, the following relation defines the pre-osteoblast to osteoblast flux:
The following relation defines the osteoblast to osteocyte flux:
The following relation defines the pre-osteoblast dynamics:
{dot over (ρ)}B*=μB*−fB*→B Eq. (3) (pre-osteoblast dynamics)
The following relation defines the osteoblast dynamics:
The following relation defines the osteocyte dynamics:
{dot over (ρ)}Y=fB→Y−ηYρY Eq. (5) (osteocyte dynamics)
Osteoclasts are multi-nucleated cells that contain many mitochondria and lysosomes. Osteoclasts cause bone resorption by secreting enzymes, such as collagenase, to release minerals, such as calcium, magnesium, and phosphate into the extracellular fluid.
In some embodiments, the following dynamic mathematical relations simulate the dynamics of the osteoclast lineage by taking into account that (a) the pre-osteoclast to osteoclast differentiation is upregulated by RANKL bound to (osteoclastic) RANK and (b) osteoclast apoptosis is upregulated by estrogen and the resorption signal.
In an embodiment of an in-silico model according to the present teachings, the following relation defines the pre-osteoclast to osteoclast flux:
The following relation defines the pre-osteoclast dynamics:
{dot over (ρ)}C*=μC*−fC*→C Eq. (7) (pre-osteoclast dynamics)
The following relation defines the osteoclast dynamics:
The systems and methods according to the present teachings also incorporate the role of sclerostin in an in-silico model of hone mineralization. Sclerostin is a protein that is encoded by the SOST gene. More specifically, sclerostin is a secreted glycoprotein with a C-terminal cysteine knot-like (CTCK) domain and sequence similarity to the DAN (differential screening-selected gene aberrative in neuroblastoma) family of bone morphogenetic protein (BMP) antagonists.
Sclerostin is produced by osteocytes and undergoes independent degradation. Sclerostin production intensity by osteocytes depends on time. Further, sclerostin production is downregulated by estrogen. The following dynamic mathematical relation simulates the estrogen dynamics by incorporating the above information.
In some embodiments, the present teachings take into account the concentration of sclerostin in the simulation of bone mineralization. Sclerostin is produced by osteocytes and its main function is to inhibit bone formation. In particular, sclerostin exerts its effects mainly by interfering with a process called Writ signaling, which plays an important role in the regulation of bone formation.
Sclerostin is produced by osteocytes and undergoes independent degradation. The intensity of sclerostin production by osteocytes depends on time. Further, sclerostin production is downregulated by estrogen. In some embodiments, the following relation is employed to define the rate of change of the sclerostin density as a function of time:
Another component of an in-silico model according to the present teachings for simulating bone mineralization is bone resorption. Bone resorption is a process by which osteoclasts degrade the tissue in bones and release the minerals and signaling factors, such as TGFβ, BMP2, PDGF and IGFs, in the bone matrix. The resorption signal is produced by osteoclasts. In some embodiments of the present teachings, it is assumed that the time scale associated with synthesis and degradation of osteoclasts is much faster than the time scale of osteoclast birth and death so that the following relation holds: r∝ρC.
More specifically, since the resorption signal acts as a regulator of bone formation and is resealed by threshold concentrations, the proportionality in the above relation can be replaced with equality. In other words, the following relation can be employed: r=ρC.
As noted above and as discussed in more detail below, the present teachings employ the dynamic mathematical relations disclosed herein to simulate bone remodeling. The gain in total bone density is proportional to the osteoblast density, and is downregulated by sclerostin and is upregulated by the resorption signal. Further, the loss in total bone density is proportional to the osteoclast density. In some embodiments of the present teachings, the following dynamic mathematical relations for bone density gain rate (Eq. (10)), the bone density loss (Eq. (11)), the total bone density (Eq. (12)), the bone mineral content (Eq. (13)), and the bone mineral density (Eq. (14)) are employed:
Although in some embodiments, an exponential dynamics of the form
It is also noted that (ϵ1) and (ϵ2) are used as dummy parameters that are introduced to provide an interface for model extensions that modulate the mineralization process.
A model according to an embodiment of the present teachings for simulating the bone mineralization process can be represented mathematically in the following form, which highlights the model's generalizable structure.
In the following S=(s, e, r) is the vector of signaling factor/hormone concentrations and M=(θromo, θromo, . . . ) represents a vector ofcompounds that describe time-dependent and time-independent model-external agents such as medications or other forms of treatment, comorbidities, lifestyle-related aspects such as physical activity, diet, et . . . .
The dynamics of precursor populations is given by the following relations:
{dot over (ρ)}C*={tilde over (μ)}C*(s)−fC*→C−{tilde over (η)}C*(s)ρC* Eq. (15) (pre-osteoclasts)
{dot over (ρ)}B*={tilde over (μ)}B*(S)−fB*→B−{tilde over (η)}B*(S)ρB* Eq. (16) (pre-osteoblasts)
{dot over (ρ)}C=αC*fC*→C−{tilde over (η)}C(S)ρC Eq. (17) (osteoclasts)
{dot over (ρ)}B=αB*fB*→B−{tilde over (η)}B(S)ρB−fB→Y Eq. (18) (osteoblasts)
{dot over (ρ)}Y=fB→Y−{tilde over (η)}Y(S)ρY Eq. (19) (osteocytes)
{dot over (s)}={tilde over (μ)}
s(S)ρY−{tilde over (κ)}s(S)s Eq. (20) (sclerostin)
{dot over (Ω)}=b+−b− Eq. (21) (bone density)
b
+={tilde over (λ)}B(S)ρB Eq. (22) (bone density)
b
−={tilde over (λ)}C(S)ρC Eq. (23) (bone density)
ϕ={tilde over (μ)}ϕ(S)−{tilde over (κ)}ϕ(S)ϕ Eq. (24) (bone mineral content)
{dot over (Ω)}m=ϕΩ Eq. (25) (bone mineral density)
where differentiation rates are of the generic form:
f
X→Y={tilde over (ω)}X(S)ρX Eq. (26) (pre-osteoclasts)
For the above incarnation of the above model, the following can hold:
One advantage of the methods and systems according to the present teachings is that the in-silico models of the present teachings can be readily extended to take into account the effect of a disease and one or more medications on the bone mineralization. By way of example, such extensions can be implemented by keeping the core structure of the model intact. In the following embodiments, an “extension” is defined as a model modification whereby
One advantage of the above procedure is that because of (ii) condition, such extensions can be modified without changing the structure of the core model. Further, because of (ii-a) condition, extensions do not change the interpretation of parameters of the core model.
In some embodiments, one or more interfaces can be provided for simulating the effect of a physiological condition (e.g., age, menopause, etc.) a disease and/or a therapeutic treatment on bone remodeling. As discussed in more detail below, in some such embodiments, the therapeutic treatment includes the administration of two or more medications to a patient. Examples of such interfaces are discussed below.
(i) Post-Menopausal Decline of Estrogen Production
It is known that estrogen levels decrease in post-menopausal women. In some embodiments, age-dependent estrogen production in women is simulated and such simulation is employed as an interface to the core model for simulating bone remodeling. More specifically, assuming fast kinetics compared to the change in production rates, the relative estrogen supply (quasi steady-state) can be given by the following relation:
where, aϵ denotes the age at the onset of estrogen decline and τe is a characteristic time scale.
In the following discussion, to compare the model output to clinical data, the model observables are related to established biomarkers, most notably markers of bone turnover.
Clinical data suggests that the biomarker serum levels change relative to the baseline scales nonlinearly with respect to formation/resorption rate. For instance, serum levels of the bone formation marker PINP (Procollagen type 1 amino-terminal propeptide) peak at about two-fold base levels 3 months after the start of sclerostin mAb treatment but BMD (bone mineral density) changes suggest an approximate 20-fold increase in the bone formation rate.
Hence, in some embodiments, a functional relationship g(b) that relates biomarker concentration M and bone formation/resorption rate τ is assumed to satisfy g (αb0)=βg(b0), where b0 represents the basal rate and the scaling factors α and β describe the change in rate and biomarker concentration after start of treatment, respectively. In particular, in some embodiments, the following relation is employed: g(b)=bq, such that b=(ln β)/(ln(α)).
In some other embodiments, a more flexible approach, which uses two parameters per biomarker, employs the following relation:
φ(b)=c+γbq Eq. (35)
where c represents basal biomarker levels and γ paramterizes the dependence of the biomarker on the respective turnover rate.
By way of example, in some embodiments, the following bone turnover markers can be employed:
(i) Bone Formation:
(ii) Bone resorption:
In this embodiment, the following relations are employed for the following biomarkers: BSAP, PINP, and CTX
φBSAP(b+)=(b+)qBSAP Eq. (36)
φPINP(b+)=(b+)qPINP Eq. (37)
φCTX(b−)=(b−)qCTX Eq. (38)
The above relations indicate that the concentrations of biomarkers BSAP and PINP depend on the bone formation rate b+, whereas concentration of the biomarker CTX depends on the bone resorption rate b−. Further, the relations indicate that three relations include different exponents.
In some embodiments, an ire-silica model according to the present teachings can provide modules for determining the effects of one or more medications on the progression of osteoporosis. In particular, in some such embodiments, such modules can be configured to determine the effect of the virtual administration of a combination of two or more medications on the progression of osteoporosis. In particular, an in-silico model according to the present teachings can take into account the cross reactions of two or more medications via capturing their mechanisms of action and determine their joint downstream effects relevant for bone turnover.
In some embodiments, the presence of a medication in a patient's body can be represented by an effective (dimensionless) variable 0 that indicates the relative concentration of a medication. Typically, the variable 0 is given in multiples of some threshold that parametrizes the effect of the drug. In general, the pharmacokinetics of a particular medication can inform the interpretation of 0.
In many embodiments, the pharmacokinetics of any drug X can be described by the following two parameters: the efficacy EX and the medication's half-life TX. A variety of different pharmacokinetic models can be employed. For example, for repeated administrations with doses c1, . . . , cn at times t1, . . . , tn, the relative concentration variable for drug X can be given by the following relation:
θX(t)=EXΣi=1nci2−(t-t
where Θ is the Heaviside step function. It should be understood that other pharmacokinetic models can also be utilized in the practice of the present teachings. For example, the above relation (39) can be replaced with other relations indicative of different pharmacokinetic models.
By way of example, the present teachings can be employed to determine the effect of RANKL inhibitors, Sclerostin antibodies, and bisphosphonates on the progression of osteoporosis as discussed in more detail below. It should, however, be understood that the present teachings are not limited to these medications, but rather can be employed to model the effects of other medications, as well as combinations of two or more medications, on bone remodeling, and in particular the progression of osteoporosis.
A RANKL inhibitor can bind the RANKL and block its interactions with RANK, which, in turn, affects pre-osteoclast differentiation.
Denosumab, a monoclonal antibody, is an example of a RANKL inhibitor that binds with affinity to RANKL. Further, without being limited to a particular theory, Denosumab appears to facilitate the secondary mineralization, which leads not only to a halt in BMD decline but also to a long-term increase in BMD. Without being limited to a particular theory, Denosumab suppresses osteocytic dissolution action, which can allow for progression of the secondary mineralization process and an increase in overall BMD.
In some embodiments, the extension of an in-silico model according to the present teachings for taking into account the effect of virtual administration of a RANKL inhibitor can be achieved by providing an interface to the variables indicating the differentiation rate of pre-osteoclasts to osteoclasts and mineralization. More specifically, these variables can be adjusted as shown below by the use of additional regulatory elements:
ωC*→ωC*−ωr−inhC*Γ+(θr−inh) Eq. (40)
ϵ→ϵ−υr−inhφΓ+(σr−inh) Eq. (41)
where θr−inh is the functional RANKL inhibitor concentration in multiples of EC50, which, for simplicity, in this embodiment, it is assumed to be identical for the regulation of differentiation and mineralization. Furthermore, ωr−inhC*<ωC* and νr−inhφ parametrize the respective maximum effect strength.
In some embodiments, the present teachings provide an interface to one or more variables that allow simulating the effect of the administration of sclerostin antibody to an individual. More specifically, in some embodiments, the effects of the administration of a sclerostin antibody can be taken into account by modifying the sclerostin dynamics by an additional regulatory element via the replacement {dot over (s)} of {dot over (s)}* and variables as follows:
{dot over (s)}→{dot over (s)}−κ
mAb
sθmAbs+δss* Eq. (42)
{dot over (s)}*=κ
mAb
sθmAbs−κ*ss*−δss* Eq. (43)
Where θmAb is the functional antibody concentration and κs is the sclerostin degradation rate; this parametrization implies that θmAb is given in multiples of the antibody concentration needed to achieve the same binding rate as the natural degradation process of sclerostin. Furthermore, s* is the concentration of antibody-bound sclerostin with ds being the unbinding rate.
To provide a single interface for multiple sclerostin antibodies, one can define the effective Romosozumab and Blosozumab doses θmAb=EX{umlaut over (θ)}X, where θX is the absolute antibody dose and the efficacy EX rescales {umlaut over (θ)}X to a dimensionless number that represents its relative contribution to the total effects of the antibodies. Thus,
θmAb=θromo+θbloso Eq. (44)
Bisphosphonates (e.g., Alendronate, Risedronate, Ibandronate, Zoledronate) inhibit osteoclast-mediated bone resorption by promoting osteoclast apoptosis.
In some embodiments, the effect of the administration of bisphosphonates on bone mineralization, and in particular, on the progression of osteoporosis can be modeled by providing an interface for adjusting the osteoclast apoptosis rate as follows:
{dot over (ρ)}C→{dot over (ρ)}C−ηbisphCΓ+(θbisph)ρC Eq. (43)
where ηbisphC is the reference additional apoptosis rate due to the presence of bisphosphonates.
In some embodiments, in order to provide a single interface for multiple bisphosphonates, effective doses θX=EX{umlaut over (θ)}X for X=alendronate, ibandronate, zoledronate, . . . , where {umlaut over (θ)}X is the absolute bisphosphonate dose and the efficacy EX rescales {umlaut over (θ)}X to a dimensionless number that represents its relative contribution to the total effects of antibodies. Thus,
θbisph=θalendro+θibandro+θzoledro+ . . . Eq. (46)
In some embodiments, a model according to the present teachings can be configured to simulate the effect of a disease on bone mineralization, and particularly, on the progression of osteoporosis. For example, in some embodiments, the present teachings can provide an interface to one or more variables that can be adjusted to account for the effects of a disease on bone mineralization.
Some examples of diseases for which the present teachings can provide an interface to the core in-silico model for simulating their effects on bone mineralization are discussed below. It should, however, be understood that the present teachings can be applied to simulate the effects of other diseases on bone mineralization as well, e.g., by providing an interface to one or more dynamic variable of the core model.
Multiple myeloma (MM) is a B-cell malignancy characterized by the clonal proliferation of malignant plasma cells in the bone marrow microenvironment. In MM, bone resorption is stimulated and new bone formation is inhibited. In particular, MM cell can stimulate osteoclast differentiation via production of IL-3, DcR3, CCL3, MIP-1β, VEGF, TNF-α, and RANKL. Further, MM cells can also adhere to bone marrow stromal cells via very late antigen 4 (VLA-4) and vascular cell adhesion molecule 1 (VCAM-1), which can in turn promote osteoclast differentiation and activation.
In this embodiment, the pre-osteoclast to osteoclast differentiation can be employed as an interface to the variable indicative of the pre-osteoclast to osteoblast differentiation rate (ωC) in the following way:
where tmm is the time of onset, τmm is the duration of the build-up and {tilde over (Θ)} is a smooth step function.
Glucocorticoids are corticosteroids, which are a class of hormones that are produced by the adrenal gland. Glucocorticoids have a role in regulating the glucose metabolism but can also modulate the immune system, where they can suppress the immune function.
Glucocorticoids can be used to treat over-functions of the immune system, such as allergies, asthma, autoimmune diseases, and after transplantation.
Glucocorticoid-induced osteoporosis (GIOP) is known to be the most common cause of secondary osteoporosis. High levels of glucocorticoids can dramatically decrease bone formation rate, osteoblast numbers and osteocyte numbers and activity.
Glucocorticoids can also increase the production of macrophage colony stimulating factor (M-CSF) and RANKL, and decrease the production of osteoprotegerin (OPG) by osteoblastic cells and osteocytes, thus resulting in an increase in both the number and activity of osteoclasts.
Glucocorticoid-induced osteoporosis can result in a rapid bone loss due to a combination of increased bone turnover and an inhibition of bone formation.
In some embodiments, the present teachings employ the reference gain rates and the differentiation rates associated with osteoblasts, osteocytes, and pre-osteoclasts as interfaces for simulating the effects of glucocorticoids on bone mineralization, and in particular their effects or the progression of osteoporosis. For example, in this embodiment, these gain and differentiation rates can be adjusted as indicated below to simulate the effects of glucocorticoids on osteoporosis:
As discussed below, in some embodiments, an in-silico model according to the present teachings can be used to implement a “Virtual Clinic.” For example, an in-silico model according to the present teachings can be used to predict the progression of osteoporosis based on data gathered from an individual patient or a collection of patients. In such an implementation, avatars (herein also referred to “digital twins”) are generated as representations of individual patients or patient populations through a set of patient-specific parameters that characterize their physiological response within the scope of a physiological mathematical model according to the present teachings.
By way of example, avatars can be used to predict patient-individual response and to optimize treatments on a personalized basis. In some embodiment, a large number of avatars, e.g., in a range of 1 to a few thousands, for example, more than 10,000, can be used to optimize general population level treatment schemes in a “Virtual Trial” setting so as to find a treatment (or a set of treatments) that collectively maximize the desired clinical outcome.
By way of example, in some embodiments, within an in-silica model according to the present teachings, avatars can be defined as follows. In particular, in many embodiments, each parameter associated with the model and its extensions can fall into the following two categories:
(1.) Parameters that are approximately constant among a patient population of interest,
(2) Parameters that reflect patient-specific physiological responses and properties.
An avatar (which is also herein referred to as a digital twin) is defined as a specific set of values for parameters in the second category. Such parameters may include, for example, the osteoporosis-related physiology of an individual patient or a patient population with similar characteristics.
For a specific patient population, an avatar may be found by fitting the model to clinical datasets (e.g., to hone mineral density measurements or measurements of biomarkers) and/or by inferring parameters from model-external properties (such as age, gender and other demographic properties). By way of example, the following datasets may be used to create avatars:
In some embodiments, the model can be based on data from individual patients. Such an approach would typically require a long history of measurements (e.g., many years, or even decades) for each individual patient. In some cases, it may be possible to infer patient-specific parameters with high confidence from auxiliary data.
In some embodiments, an in-silico model according the present teachings allows simulating the effect of a physical activity and/or altered gravitational conditions in bone remodeling, and in particular in treatment and/or progression of osteoporosis.
For example, in a model according to an embodiment, the effects of a time-dependent variation in physical activity can be incorporated by a dynamical inhibition of sclerostin. Let ϑx (t) denote a coarse-grained normalized physical activity, e.g., fraction of a month spent on exercises of type x (e.g., running, cycling, weight lifting), such that ϑx=1 corresponds to the average of a considered patient population. A corresponding module modulates the sclerostin production according to
μs→μs+2μphyss(Γ−(ϑx)−½) Eq.(52)
where μphyss<μ
Similarly, the effects of altered gravitational conditions, e.g., higher or lower gravity, can be incorporated into an in-silico model according to the present teachings. For example, a relation similar to the above Eq. (52) can be employed in which another variable ϑx(t) with type x indicating gravity as an effect can be employed, so that ϑx=1 corresponds to normal conditions (i.e., gravity on the earth's surface), ϑx<1 corresponds to lower gravity and ϑx>1 corresponds to higher gravity.
In some embodiments, the present teachings allow the selection of an optimal protocol for treatment of osteoporosis (e.g., slowing the progression rate of osteoporosis). By way of example, such selection of an optimal protocol can include the selection of an optimal dosage and/or administration regimen of a drug. With reference to
As discussed below, these modules are in communication with one another to provide in-silico simulation of bone remodeling, and in particular, progression of osteoporosis using the dynamic mathematical relations discussed above. In particular,
With reference to both
By way of example, and without limitation, a cost function may be defined as a weighted sum of residuals. For example, let Ωm,i and
D
1(P)=Σi=1n(Ωm,i−
D
2(P)=Σi=1n(φiP1NP=
D
3(P)=Σi=1n(φiCTX=
Then, one possible cost function is given by
C(P)=Σκ=13wκDκ(P) Eq. (56)
where the wi denote the weights of the respective quantity in the cost function.
An alternative cost function can be defined as a sum of logarithms, which renders the cost function less dependent on the absolute scale of respective quantities,
C(P)=Σ=κκ=13wk ln Dκ(P) Eq. (57)
Other cost functions can also be employed in the practice of the present teachings.
The fitting module is configured to fit the data to a predefined model so as to derive parameters that will be employed by the avatar module 1004. Some examples of model parameters include, without limitation, thresholds for regulation of bone turnover through estrogen, sclerostin; rates of cell birth, differentiation and death; rates of bone formation and resorption per cell density. Typical avatar parameters can include efficacy and half-life of the administered drugs, parameters determining the relation between bone turnover and bone turnover markers found in a patient, time scale of estrogen decline after onset of menopause, etc. The fitting module then formats the data so as to enable its comparison with simulations and sets up a parameter fit strategy.
To fit the parameters, the fitting module 1006 invokes the avatar module to compute thousands of avatar realizations and compares them with the data. For each realization, the fitting module 1006 passes a candidate parameter set (i.e., avatar and model parameters) to the avatar module and the avatar module computes the entire lifetime of bone evolution for that avatar.
For example, as shown in the workflow chart of
Subsequently, the avatar module 1004 passes this information, together with the entire set of parameters, to the simulation module 1002, which employs a numerical solver to compute a solution, of the bone remodeling (osteoporosis) model according to the present teachings using all parameters and auxiliary functions obtained from the avatar module 1004. In particular, as depicted in the workflow of
The simulation module 1002 then passes this information back to the avatar module 1004. The avatar module 1004 employs this information to compute additional avatar-specific quantities, such as time-dependence of the bone turnover markers. In general, in many embodiments, the avatar module handles all the characteristic that are specific to the patient.
The avatar module 1004 then compiles all the results and passes the compiled information back to the fitting module 1006. The fitting module 1006 in turn computes the value of a cost function, which had been passed earlier to the fitting module as discussed above, between the avatar results and the clinical data. If the calculated cost function, or a reduction in calculated value of the cost function, is below an acceptable threshold, the avatar calibration is completed. However, if the calculated cost function, or a reduction in calculated value of the cost function, is not below acceptable threshold, the fitting module 1006 will generate new parameter set based on the previous results, and passes the new parameters to the avatar module. This process is repeated iteratively until an acceptable convergence of the model to the clinical data is achieved. In some cases, it is possible for the value of the best (i.e., the optimal) cost function to be greater than a preset threshold. In such a situation, the avatar associated with that cost function can be rejected. In this embodiment, the fit parameters vary between different candidate sets while all other parameters remain fixed to a reference value.
The workflow for obtaining an optimal virtual treatment using an osteoporosis model according to the present teachings can vary depending on the type of scenario under investigation. For example, in some cases, it may be desired to optimize a dosing scheme for a single drug. In such a case, there are different candidate dosing schemes (e.g., administering the same total drug dose but changing the dosage and the intervals between the administration of the drug). The system 1000 can be used to run each scheme for a population-level avatar or for an entire set of different patient-level avatars (i.e., a so-called Virtual Trial). Population-level avatars can be generated, for example, based on averages of parameters associated with members of a population. The probed dosing schemes can be ranked in accordance with their performance. For instance, in some cases, the desired goal may be the maximization of gain in bone mineral density, or the desired goal may be to optimize the temporal period after a treatment during which the bone density gain due to the treatment is sustained, or combinations thereof.
In some cases, the system 1000 can be employed to probe combination therapies. The workflow is similar to that discussed above, where certain dosing schemes (e.g., sequential and parallel administration of different drugs in different relative doses/intervals, etc.) can be probed using the generated avatars and ranked according to a pre-defined performance criterion. By way of example, such criteria may be based on maximizing the time until BMD falls below the baseline level at the start of the treatment, and/or minimizing potential side effects of a drug. In some embodiments, the criteria can be based on maximizing the total BMD gain (e.g., maximum BMD at some point during the treatment). In some embodiments, each of the simulation module 1002, the avatar module 1004, and the fitting module 1006 can be implemented in software, where each module is implemented as a software layer. For example, in some implementation, an object-oriented programming language, such as Java® or Python®, can be employed and each of the simulation module 1002, the avatar module 1004, and the fitting module 1006 can be implemented as a class, which can be instantiated as an object during execution.
The software modules can be executed on any suitable hardware platform. By way of example,
By way of example, the processor 400 can be a general and/or special purpose microprocessor, such as an application-specific instruction set processor, graphics processing unit, physics processing unit, digital signal processor, image processor, coprocessor, floating-point processor, network processor, and/or any other suitable processor that can be used in a digital computing circuitry. Alternatively or additionally, the processor can comprise at least one a multi-core processor and a front-end processor. By way of example, the permanent memory 402 can be, for example, magnetic disks (e.g., internal or removable disks), magneto-optical disks, one or more of a semiconductor memory device (e.g., EPROM or EEPROM), flash memory, CD-ROM, and/or DVD-ROM disks.
Instructions and data for performing in-silico simulations according to the present teachings can be stored in the permanent memory and can be transferred onto the RAM during execution.
The following examples are provided for further elucidation of various aspects of the present teachings and are not intended to necessarily indicate the optimal ways of practicing the present teachings and/or the optimal results that may be obtained.
An example of an implementation of methods and systems according to the present teachings for in-silico simulation of osteoporosis was employed to obtain the simulation results depicted in
Parameters of the osteoporosis model were obtained using a simultaneous numerical fit of the model including modules for sclerostin antibodies and RANKL inhibitors to clinical datasets. The following clinical datasets were employed for obtaining the depicted simulation results: Recknor et al., J. Bone Miner. Res. 30 (9), 1717-1725 (2015) (
A fit of the model to the clinical data performed in a manner discussed above resulted in the following parameter list:
Those having ordinary skill in the art will appreciate that various changes can be made to the above embodiments without departing from the scope of the invention.
The transient population burst through a step-like input can be computed as discussed below.
Consider a progenitor population x and an active/mature population y, subject to the dynamics
{dot over (x)}=μ−ωϕ(t)x
{dot over (y)}=ωϕ(t)x=κy
where ϕ is an external input. Initially, ϕ=1 and x and y are at steady state, x0=μ/ω and y0=μ/κ. At time t=0, ϕ jumps to a new value ϕ0>1, so that
i.e., x evolves from the old steady state x0 to the new steady state x0/ϕ0. Hence, y evolves as
i.e., as two competing exponential decays that produce a transient increase in y, followed by a return to y0. We use this mechanism as the explanation for the transient effect of Sclerostin antibodies on bone formation, where the role of x is played by the pre-osteoblasts and the role of y is played by the osteoblasts.
