Model checking the biological model of membrane computing with probabilistic symbolic model checker by using two biological systems

Problem statement: Membrane computing formalism has provided better modeling capabilities for biological systems in comparison to conventional mathematical models. Model checking could be used to reason about the biological system in detail and with precision by verifying formally whether membrane computing model meets the properties of the system. Approach: This study was carried to investigate the preservation of properties of two biological systems that had been modeled and simulated in membrane computing by a method of model checking using PRISM. The two biological systems were prey-predator population and signal processing in the legend-receptor networks of protein TGF-β. Results: The model checking of membrane computing model of the biological systems with five different properties showed that the properties of the biological systems could be preserved in the membrane computing model. Conclusion: Membrane computing model not only provides a better approach in representing and simulating a biological system but also able to sustain the basic properties of the system.


INTRODUCTION
Membrane computing provides a hierarchical structure for molecular computation in which embraces play an essential role for objects to pass in a regulated fashion within and across the membranes. Since its inception in 1998, much research on theoretical aspects has been done to establish membrane computing computational power (Paun, 1998;2000). In recent time, the research interest is more concentrated in applying the membrane computing formalism to solve real world problems. One of the attempts is using membrane computing capabilities in modeling biological systems. Studies of biological systems such as cells in silico have greatly reduced the need for expensive and prolonged lab experiments. Construction of a biological system into the membrane computing model provides a better understanding of the dynamics and functionality of the system. The biological description of membrane computing formalism has been utilized to characterize and preserve the elements in biological systems. The research in this line shows that biological system can be modeled better using membrane computing than the conventional methods using mathematical representation such as Ordinary Differential Equation (ODE) (Bernardini et al., 2006;Muniyandi and Abdullah, 2009;2010). Simulation strategies based on membrane computing such as Gillespie Algorithm are also used to proof this claim.
Meanwhile, through simulations the operation of the model can be studied and consequently, the properties concerning the behavior of the biological system can be inferred as shown by Muniyandi and Abdullah (2010). However, although simulation outlines some key features, the preservations of other properties of model such as computability, understandability, extensibility and relevancy also need to be checked to make sure essential properties, dynamic behaviors, informal concepts and higher levels extensible capabilities of specific biological system are captured as well (Regev and Shapiro, 2004).
In comparison to deterministic models, stochastic models are more difficult to visualize and due to that it is necessary to reason about the systems in detail and with precision. Model checking is a verification procedure to verify formally whether a model meets the properties stated in the specification. The system's specification is represented by using temporal logic formulas and symbolic algorithms are used to navigate the model to check if the specification is fulfilled or not.
In recent time with the rapid research developments in Systems Biology, the use of model checking for the analysis of biological models has attracted much attention. Some of the research are: The analysis of transcriptional regulation by LTL model checking (Barnat et al., 2009); Performing statistical model checking using Bayesian sequential hypothesis Testing on biological systems (Sumit et al., 2009); algorithmic algebraic model checking on metabolic networks (Mysore and Mishra, 2007); temporal logic analysis of Gene Networks (Batt et al., 2007) and probabilistic model checking of complex biological pathways (Heath et al., 2006;Kwiatkowska et al., 2008).
There are also studies investigating the use of model checking for membrane computing models such as natural algebraic specification for the membrane computing (Andrei et al., 2005) investigate various models of membrane computing to identify their model-checking problems (Dang et al., 2005) and attempt to use PRISM (probabilistic and symbolic model checker) with a membrane computing model .
In this study, we investigate whether the properties of two biological case studies modeled using membrane computing and simulated by Muniyandi and Abdullah (2010) are preserved by using the PRISM model checker. The two biological case studies are preypredator population and signal processing in the legendreceptor Networks of protein TGF-β. These two case studies have been modeled in membrane computing formalism and simulated with simulation strategy of membrane computing with Gillespie Algorithm (Muniyandi and Abdullah, 2009;2010).

MATERIALS AND METHODS
PRISM model checker: PRISM (Kwiatkowska et al., 2002) is probabilistic model checker that represent a technique for formally verifying quantitative properties of a stochastic system. A model to be analyzed is specified in PRISM language, which is a simple, high level, state based language based on the Reactive Modules formalism (Alur and Henzinger, 1999). PRISM supports a probabilistic model based on Continuous Time Markov Chains (CTMC) and systems specifications through continuous stochastic logic for stochastic systems.
Case studies: Prey-predator population: The primary example of a Prey-Predator population (Jones and Sleeman, 2003) is a system comprised of a plant population and an herbivorous animal dependent on that plant for food. In this case, the predator is totally dependent on the prey as its only food supply. The prey species has an unlimited food supply and no threat to its growth other than the specific predator. If there were no predators, the prey species grows exponentially. But there are predators, which must account for a negative component in the prey growth rate. The assumptions for the model are the rate at which predators encounter prey is jointly proportional to the sizes of the two populations and a fixed proportion of encounters lead to the death of the prey. Based on these behaviors, three rules are formulated in membrane computing as follows (k 1 , k 2 and k 3 are kinetic constants):

Properties: The properties for the Prey-Predator
Population and TGF-β are obtained from the behavior of these systems elaborated in the respective research studies (Jones and Sleeman, 2003;Villar et al., 2006).

Properties of prey-predator population:
Recurrence behavior of the system and the existence of equilibrium probability distribution maintain the stability of the system in the form of oscillations. To facilitate this behavior, the prey-predator system should preserve the following properties: • The rules are selected stochastically based on the number of preys and predators at each time steps and the value of reaction constants to maintain the equilibrium of the system • The number of preys and predators must not equal to 0 at any time steps • The number of preys and predators become equal or intersect each other twice at each cycle of the system • Percentage of increase or decrease of number of preys are higher most of the time steps than the number of predators • Percentage of change between preys and predators are higher most of the time steps for preys than predators Properties of TGF-β β β β: The purpose of this model is to study the signal processing potential of the ligandreceptor network and receptor trafficking. The signaling activity is proportional to the number of legendreceptor complexes that are present in the internalized endosomes. In order to attain this behavior, the other essential properties required for this model are: • Legends induce the formation of receptor complexes with type I and type II receptors • Receptors and legend-receptor complexes can be present in two spatially distinct compartments: Plasma membrane and internalized endosomes • Receptors and legend-receptor complexes are continuously internalized into endosomes and recycled back to the plasma membrane • Receptor degradation has a constitutive contribution, which is the same for free receptors and legend-receptor complexes • Receptor degradation has a legend-induced contribution, which affects only receptors that have been complexes with legends The biological case studies modeled in membrane computing by Muniyandi and Abdullah (2010) are translated into PRISM formalism. This translation technique from membrane computing into PRISM applied in this research is proposed by . Then the model in PRISM is simulated and model checked with the properties for each of the case studies. PRISM is used to specify and to analyze properties based on rewards. Rewards are used to reason the behavior of the model in a certain fashion by measuring the probability as well as to identify a wide range of quantitative measures relating to modeling behavior. The data of the results generated by PRISM is exported into Scilab (Campbell et al., 2000) to be presented into a graph. Finally, the graph is analyzed to verify whether the identified properties are preserved or not. endmodule By using the concept of rewards, PRISM is used to specify and to analyze properties of Prey-Predator as follows:

Model
Property 1: Rewards "reactions" is to analyze the stochastic behavior of the system to maintain the equilibrium. The rewarded value of 2, 4 and 6 refer to the selection of rules R1, R2 and R3 respectively at each time steps: rewards "reactions" [R1] true: 2; [R2] true: 4; [R3] true: 6; endrewards The Fig. 1 shows four graphs indicating the selection of rules at four different periods of time steps based on the rewards. These graphs show that at each time step one of the three rules is selected stochastically. The patterns of selections differ for each period of time steps as shown by each graph. This means that the stochastic behavior of the system is maintained to make sure the stability and consistency of the system at each cycle of the oscillation. Property 2: Rewards "PreyPredatorgreater0" is to verify that the number of preys or predators is always greater than 0 throughout the period of oscillation. The rewarded value of 1, 2 and 3 are assigned to each of the following conditions respectively, when the condition is true at each time steps: Both prey and predator equal to 0; prey greater than 0 and predator equal to 0 and prey equal to 0 and predator greater than 0. If the three conditions are not met, then the rewarded value is 0: rewards "PreyPredatorgreater0 " Prey = 0 & Predator = 0: 1; Prey>0 & Predator = 0: 2; Prey = 0 & Predator>0: 3; endrewards Figure 2 shows that not at once along the simulation steps of the system, the number of prey or predator has been equal to 0 based on the rewards. This result demonstrates that the behavior of the system is always consistent to make sure that the number of prey and predator must always greater than 0 to stabilize the system. Property 3: Rewards "Intersection" is to determine the expected number of intersection at each cycle of the oscillation. Rewarded value of 1 is generated when the number of preys and predators are equal that is there is an intersection between prey and predator: rewards "Intersection" Prey = Predator: 1; endrewards  Figure 3 demonstrates the time steps when number of prey is equal to number of predator generated by the rewards. The intersection between prey and predator occurs twice at each cycle in the oscillations when the prey and predator either keep decreasing or increasing. However the graph shows that the period of occurrence of one intersection to another is not similar. This is mainly due to the stochastic behavior of the system. However the pattern of intersection in each of the cycle in the oscillation is preserved to maintain the stability of the system. Property 4: Rewards "PreyIncreaseDecrease" and "PredatorIncreaseDecrease" are to measure the percentage of increase or decrease of prey and predator from their respective initial amount, at each time step: rewards "PreyIncreaseDecrease" true: ((Prey-1000)/1000)*100; endrewards rewards "PredatorIncreaseDecrease" true: ((Predator-200)/200)*100; endrewards Figure 4 illustrates the percentage of increase and decrease of prey and predator based on the rewards. The graph shows that the percentage of increase and decrease of predator is higher than the percentage of increase and decrease of prey. At the initial state, the number of prey is five times higher than the number of predator. The sharp increase of predator is to certain extent decrease the population of prey. Meanwhile, the sharp decrease of predator population gives some space for prey to increase its population to attain the initial level over again. The percentage of increase/decrease of prey and predator is almost similar at each of the cycle of the oscillation. This demonstrates that the equilibrium of prey-predator population has been preserved by maintaining the percentage of increase/decrease of prey and predator accordingly at each time step.

Fig. 5: Percentage of changes between preys and predators
Property 5: Rewards "PreyChanges" and "PredatorChanges" are to measure the percentage of differences between the amount of preys and predators in the population, at each time step: rewards "PreyChanges" true: ((Prey-Predator)/Prey)*100; endrewards rewards "PredatorChanges" true: ((Predator-Prey)/Predator)*100; endrewards Figure 5 shows the percentage of changes of prey and predator compared to the opposite population as outlined in the rewards. This graph shows that at most of the time number of prey exceeds number of predator. But at certain period when number of predator is above the number of prey, prey is decreasing to control the increase of predator. When number of predator is decreasing, the number of prey is increasing. As shown in the Fig. 4, Prey Predator these results also show that the stability of the system is preserved by controlling the increase and decrease of the prey and predator accordingly at each time step. module plasma_membrane Property 1: Rewards "CF" is to measure the process of legends inducement in the formation of receptor complexes (IRIRII1) with type I (RI) and type II (RII) receptors. This process is performed by reaction CF in the PRISM model and rewards "CF" measures the concentration of IRIRII1 as the product of this reaction activation at each time step: rewards "CF" [CF] true:IRIRII1; endrewards Figure 6 shows the execution of reaction CF based on rewards "CF". At the initial stage this reaction is activated more regularly to meet the formation of IRIRII in plasma membrane. The increasing amount of IRIRII at this stage shows that internalization into endosomes is occurring consistently but slowly. This shows by the slow increase prey predator of IRIRII in endosomes as shown in the simulation done by Muniyandi and Abdullah (2010). Internalization is rapidly performed between time steps of 2000 and 3000 in which period the IRIRII in endosome achieve the peak. After this process, the formation of IRIRII in plasma membrane is stabilized with the combination of IRIRII recycling into plasma membrane and IRIRII internalization into endosomes.

Property 2:
The following rewards are to measure receptors and ligand-receptor complexes availability in Plasma Membrane (PM) and Endosome (E): rewards "RIAvailability_PM" true: RI1; endrewards rewards "RIAvailability_E" true: RI2; endrewards rewards "RIIAvailability_PM" true: RII1; endrewards rewards "RIIAvailability_E" true: RII2; endrewards rewards "IRIRIIAvailability_PM" true: IRIRII1; endrewards rewards "IRIRIIAvailability_E" true: IRIRII2; endrewards  However there is interaction between objects in these two compartments to generate IRIRII in endosomes to measure signal processing. This can be seen when RI and RII are internalized from plasma membrane, the concentration of these receptors is keep increasing in endosomes. So is for IRIRII until it achieves the peak. But after the peak, the concentration of receptors and legend receptor complexes are stabilized to a certain level with recycling and internalization activity.

Property 3 and 4:
Rewards "IN" is to measure the internalization of receptors and legend-receptor complexes into endosomes. The reactions I1, I2 and I3 involve in this process. Rewards "RC" is to measure the recycling back of receptors and legend-receptor complexes into the plasma membrane. The reactions RC1, RC2 and RC3 involve in this process:  Figure 9 and 10 shows receptors and ligand-receptor complexes are continuously internalized into endosomes and recycled back to plasma membrane based on the rewards. This result shows that the internalization activity is slow at the initial stage before becoming intense at the period of achieving the peak for IRIRII in endosomes. Meanwhile there is hardly any recycling activity at the initial stage but this activity gain momentum after IRIRII in endosomes achieve the peak. At last both activities become constant to stabilize the system.    Figure 11 demonstrates the reactions that activate receptors degradation based on the rewards. The graph for CD1, CD3 and CD4 shows that receptor degradation has a constitutive contribution, which is the almost similar for the receptors and legend-receptor complexes. And graph CD2 shows that receptor degradation has a ligand-induced contribution that affects legend-receptor complexes.

DISCUSSION
Membrane computing has shown that it is capable of abstracting the structure and processes of biological system to represent them in a formal way without disregarding its biological characteristics. This allows the biological systems to be modeled in better way to counter some of the limitations in the conventional method based on ordinary differential equations. Moreover the membrane computing model is capable of preserving the stochastic elements of biological systems which is absent in the ODE model.
The results of the model checking of membrane computing models for two biological systems demonstrates that the properties of the system has been preserved in membrane computing model. The model checking with Prism shows that the properties of Prey-Predator Population have been preserved in the membrane computing model by maintaining the stability and consistency of the system. This means that the stochastic behavior of membrane computing has conserved the equilibrium of the oscillation of the Prey-Predator system by controlling the growth and reduction of the population of Preys and Predators accordingly. The model checking of Legend-Receptor Networks of protein TGF-β shows the formation, degradation, internalization, recycling and synthesis processes in the system have been activated according to the properties of the system in the membrane computing model. This enables the system to perform the signaling activity according to the number of ligand-receptor complexes in the internalized endosomes.

CONCLUSION
This study demonstrates that the non-determinism and stochastic behavior of membrane computing capable in preserving the properties of the biological system that had been modeled using deterministic approach of ODE. This reinforces that the membrane computing model not only provide better approach in representing and simulating biological system but also able to sustain the properties of the original model.