Application of Petri nets in bone remodeling.

Understanding a mechanism of bone remodeling is a challenging task for both life scientists and model builders, since this highly interactive and nonlinear process can seldom be grasped by simple intuition. A set of ordinary differential equations (ODEs) have been built for simulating bone formation as well as bone resorption. Although solving ODEs numerically can provide useful predictions for dynamical behaviors in a continuous time frame, an actual bone remodeling process in living tissues is driven by discrete events of molecular and cellular interactions. Thus, an event-driven tool such as Petri nets (PNs), which may dynamically and graphically mimic individual molecular collisions or cellular interactions, seems to augment the existing ODE-based systems analysis. Here, we applied PNs to expand the ODE-based approach and examined discrete, dynamical behaviors of key regulatory molecules and bone cells. PNs have been used in many engineering areas, but their application to biological systems needs to be explored. Our PN model was based on 8 ODEs that described an osteoprotegerin linked molecular pathway consisting of 4 types of bone cells. The models allowed us to conduct both qualitative and quantitative evaluations and evaluate homeostatic equilibrium states. The results support that application of PN models assists understanding of an event-driven bone remodeling mechanism using PN-specific procedures such as places, transitions, and firings.


Introduction
Bone remodeling is a homeostatic process for maintenance of healthy bone through a removal of old bone by osteoclasts followed by deposition of new bone by osteoblasts. 1 Since the process is a complex interplay among many molecules and varying types of cells, building a mathematical model is useful for developing therapeutic strategies for patients with metabolic disorders such as osteoporosis. Pioneering modeling works using a set of ordinary differential equations (ODEs) includes evaluation of the effects of parathyroid hormone (PTH) and the role of a signaling system composed of osteoprotegerin (OPG), receptor activator of nuclear factor κB (RANK), and RANK ligand (RANKL). [2][3][4] Although describing a bone remodeling process using ODEs is a basis for quantitative analyses, numerical results often present a challenge in interpreting the role of discrete molecular and cellular events involved in bone remodeling.
Here, we applied Petri nets (PNs) as an alternative approach to simulate interactive events in bone remodeling. Compared to the approach with ODEs, PNs offer several advantages. First, PNs provide graphical representation of individual interactions in the system that seems appropriate for modeling, analysis and simulation of large-scale dynamic systems. [5][6][7] Second, since system behaviors in PNs are monitored by discrete events through the firing of transitions, PNs can offer a framework for implementing complex temporal inter-related events (both synchronous and asynchronous) as well as structural interactions. Third, since any events are generated and transformed in the network model, not only deterministic but also stochastic processes can easily be built in. Fourth, PNs enable us to conduct both qualitative systems analysis (structural characterization) and quantitative analysis (monitoring dynamic behaviors).
PNs have been extensively applied in many engineering areas such as manufacturing systems, [8][9][10] transportation systems, 11,12 and communication networks. 13,14 Recently, PNs have been applied for modeling and analysis of metabolic pathways. For instance, qualitative analyses have been conducted focusing on place invariants and transitions invariants 15,16 as well as steady states of metabolic pathways. 17,18 Quantitative analyses have also been performed in calculation of the probability distribution of molecular species [19][20][21] and molecular concentrations. [22][23][24] Few studies, however, have been directed to both qualitative and quantitative analyses with reference to the ODE-based approach. Our particular interest herein is to evaluate ODE-driven equilibrium states using a PN model. This evaluation is especially important for physiological processes like bone remodeling, where variations from homeostatic equilibria may be linked to metabolic disorders.
In order to examine a potential capability of PNs in bone remodeling, it is neither feasible nor desirable to attempt to build models that include many unknown factors. We thus focused on one of the major signaling pathways (OPG-RANK-RANKL pathway) with four dominant bone cell types (two types of osteoblasts and two types of osteoclasts). Using PN models, our qualitative analysis was focused on identifying two properties (place invariants and transition invariants). Place invariants are for characterizing relationships among variables, while transition invariants are for identifying a set of sub-networks in the overall network. In the quantitative analysis, we evaluated the homeostatic equilibrium states based on PNs and compared them with the results obtained from ODEs.

Derivation of ODes
The bone remodeling process was modeled using 8 state variables (4 in the molecular level, and 4 in the cellular level) (Fig. 1). In the molecular level, 4 state variables focusing on OPG/RANK/RANKL pathway were x O (t), x L (t), x OL (t), x KL (t), which corresponded to the concentrations of OPG (O), RANKL (L), RANK (K), OPG-RANKL complex (OL), and RANK-RANKL complex (KL). The first-order ODEs were (  x = time derivative of state x; k i = rate constant; p i = synthesis rate; and d i = degradation rate): In the cellular level, 4 state variables represented the numbers of 4 different types of cells (OBP = osteoblast precursors; AOB = active osteoblasts; OCP = osteoclast precursors; and AOC = active osteoclasts). Amplification and differentiation of those cells were modeled: where N = number of cells; α i = synthesis rate; β i = interaction factor, and γ i = degradation rate. The predicted values of the above parameters, employed in this study, are listed in Table 1.

Modeling strategy with Pns
To illustrate our modeling strategy, a simplified version of PN models is shown (Fig. 2). PNs are weighted bipartite graphs with two types of nodes (places and transitions) and arcs. Places (circles in Fig. 2; e.g. molecules and cells) indicate the conditions under which transitions can occur, and transitions (bars in Fig. 2) mark events that alter states (e.g. synthesis, degradation, and chemical reaction). Arcs (arrows in Fig. 2) capture casual relationships as well as interactions among nodes, and they are associated with integer weights that regulate events. The states of PNs are defined by tokens (black dots inside places in Fig. 2), which represent the number of resources (e.g. number of molecules). In the example in Fig. 2, two molecules A and one molecule B are required to synthesize two molecules C, and this event is regulated by the firing of the transition t 1 .

Qualitative Pn analysis
In qualitative PN analysis, two behavioral properties (place invariants and transition invariants) were examined. Place invariants are a set of places where the number of tokens in those places remains constant during the evolution (dynamic behaviors) of the system. They identify the processes in which the numbers of molecules or cells stay unchanged. Transition invariants are a set of transitions where their sequences of firings can be reproduced in the specific states. They are useful to capture cyclic reaction processes and can be used to identify reversible subnetworks in metabolic networks.

Quantitative Pn analysis
PNs represent discrete state and event-driven systems, and quantitative PN analysis was conducted using numerical integration. In our analysis, the dynamic  behaviors in PNs were characterized using a flow of tokens triggered by firings of transitions. The regulatory rules were derived for each firing of transitions from the parameters in ODEs, and tokens were added or removed based on the network structures. Note that since ODEs offer continuous quantities, those continuous quantities were discretized in the simulation step (in terms of the number of tokens in places). Using a set of initial conditions, we traced the numbers of tokens in the places and evaluated their temporal alterations with reference to the equilibrium states derived from ODEs. Although the time axis was defined in terms of the event-driven firing sequences, it was uniquely linked to real time.

Pn model
The overall PN model for the selected bone remodeling process is illustrated (Fig. 3).  Prior to qualitative and quantitative analyses, we evaluated the sensitivity of the equilibrium states to the selected parameters. We first obtain the equilibrium states analytically from the set of ODEs: Then a partial derivative of all equilibrium states was derived with respect to each of the chosen parameters such as , , … OC γ . There were 144 derivatives corresponding to 8 state variables and 18 parameters, and 43 derivatives were non zero. Two parameters (p L and d L ) were involved in the equilibrium states of 6 state variables, while 5 parameters (α 2 , α 4 , β 2 , γ 2 , γ 4 ) were linked to a single equilibrium state only. Among 6 state variables affected by p L , for instance, the most sensitive state to a variation of p L value was x OL .

Qualitative analysis
The structural PN model in Fig. 3 did not have any place invariants. This result implies that no conservation of molecules or cells that were involved in this metabolic process. However, 11 transition invariants were identified (Table 2). For instance, the invariant sub-network p O → p 1 → k 1 → β 2 → γ 2 corresponds to a process of the synthesis of OPG-RANKL complexes from OPG and RANKL and its interaction with active osteoblasts, while the invariant sub-network p 1 → ak 3 → β 3 → α 4 → γ 4 corresponds to the interactions among RANKL, RANK-RANKL complex, osteoclast precursors, and active osteoclasts. Furthermore, the sub-networks k 2 → k 1 and ak 3 → k 4 corresponds to the reversible processes among OPG, RANKL and OPG-RANKL complexes, and between RANKL and RANK/RANL complexes, respectively.

Quantitative analysis
Simulation of a sub-network i We first examined the transient responses of a simplified PN model (sub-network I) (Fig. 4). In this sub-network, the concentration of OPG, x O (t), assigned in place p 1 , was expressed in ODE:

Transition invariants Remarks
Synthesis of OPG-RAnKL from OPG and RAnKL, and its interaction with active osteoblasts

Synthesis and degradation of OPG
Closed-loop interactions among RAnKL, RAnK-RAnKL, and osteoclast precursors Simulation of the sub-networks i and ii We next evaluated the interaction between two sub-networks, which were described in ODEs as: and  x t p d x t Evaluation of the equilibrium states using the entire Pn model The transient responses for the entire PN model were simulated using the initial conditions that deviated from the ODE predicted equilibrium states. Although time required for reaching steady states varied among 8 state variables, all 8 variables returned closely to the ODE equilibrium states (Figs. 6 and 7). First, the results for the molecular network in Fig. 6 revealed that the steady state PN values were 571, 3, 1633, and 89 nM for x O , x L , x OL and x KL , respectively. The ODE predicted values were 571, 2.86, 1630, and 85.7 nM in this order. Second, the steady state PN values for the cellular network in Figure 7 exhibited 79,975 for N OBP , 40 for N AOB , 9 for N OCP , and 1 for N AOC , while the ODE predictions were 80,000, 40, 9.0, and 0.9 respectively.
In the current study, we conducted both qualitative and quantitative analyses. Our qualitative analysis allowed us to identify reversible processes, and determine interactions and dependencies among molecules and cells. The quantitative analysis for equilibrium states enabled to establish a bridge between the ODE-based continuous responses and the event-driven discrete networks. The potential capability of PNs in investigating metabolic networks is multifold. First, unlike ODEs PN models can easily incorporate non-differentiable functions. For instance, administration of therapeutic agents such as OPG can be given in an arbitrary form including a series of impulsive dosages. Second, the potential effect of individual molecules such as RANK can be monitored graphically in any sub-networks. Third, an effect of a single event (e.g. synthesis of one molecule) in the entire PNs can be evaluated. Fourth, differential transient responses and time constants can be determined through temporal evolutions among variables. Lastly, although the described bone remodeling model is much simpler than a true physiological phenomenon, the present PN model can be expanded by adding more places and transitions.
Since OPG can reduce bone resorption through interactions with RANK and RNAKL, it can be used as a therapeutic agent for patients with osteoporosis. 25 In order to achieve a suitable outcome without inducing potential side effects, the administration sequence (timing and dosage) needs to be evaluated. We believe that the PN model in the current study can be used to predict a safe, effective administration strategy. Bone is a complex organ, and biological and mechanical characters differ depending on locations.
Although the current PN model does not include those local variations, it is possible to differentiate site-specific dynamics by considering additional state variables and parameters.

conclusion
The study herein presented a unique PN model for evaluation of bone remodeling focusing on the OPG/ RANK/RANKL signaling pathway among precursor/ active osteoblasts and osteoclasts. The described PN model is able to characterize qualitative (structural) and quantitative (dynamic) properties of the complex homeostatic process. It identified transition invariants (closed-loops and reversible processes) and verified the equilibrium states derived from the associated set  publish with Libertas Academica and every scientist working in your field can read your article "I would like to say that this is the most author-friendly editing process I have experienced in over 150 publications. Thank you most sincerely." "The communication between your staff and me has been terrific. Whenever progress is made with the manuscript, I receive notice. Quite honestly, I've never had such complete communication with a journal." "LA is different, and hopefully represents a kind of scientific publication machinery that removes the hurdles from free flow of scientific thought." Your paper will be: • Available to your entire community free of charge • Fairly and quickly peer reviewed • Yours! You retain copyright http://www.la-press.com of ODEs. Since PN's discrete network modeling fits well to event-driven physiological responses, further application of PN models should contribute to the understanding of complex molecular and cellular interactions and development of therapeutic strategies in bone remodeling and other biological processes.

Disclosure
The authors report no conflicts of interest.