Gexiang Zhang

Miguel Á Martínez

Mario J Pérez-Jiménez

Zeyi Shang
email: zeyi.shang@lacl.fr

Sergey Verlan
email: verlan@u-pec.fr

An Overview of Hardware Implementation of Membrane Computing Models

Keywords: CCS Concepts:, Hardware → Simulation and emulation, Transaction-level verification, Semi-formal verification Membrane Computing, P Systems, Hardware Implementation, Field Programmable Gate Array (FPGA), Compute Unified Device Architecture (CUDA), Graphic Processing Unit (GPU)

The model of membrane computing also known under the name of P systems is a bio-inspired large-scale parallel computing paradigm having a good potential for the design of massively parallel algorithms. For its implementation it is very natural to choose hardware platforms that have important inherent parallelism, like field-programmable gate arrays (FPGAs) or compute unified device architecture (CUDA)-enabled graphic processing units (GPUs). This paper performs an overview of all existing approaches of hardware implementation in the area of P systems. The quantitative and qualitative attributes of FPGA-based implementations and CUDA-enabled GPU-based simulations are compared to evaluate the two methodologies.

INTRODUCTION

With the arising of protocells 3.84 billion years ago in hydrothermal vent precipitates [START_REF] Dodd | Evidence for early life in Earth's oldest hydrothermal vent precipitates[END_REF], the evolution of unicellular organisms led to the apparition of multicellularity [START_REF] Pohorille | Self-assembly and function of primitive cell membranes[END_REF]. The biological cell structure defined by the membranes has evolved and was optimized for billions of years. As a consequence, a biological cell is a powerful parallel processing unit that can perform sophisticated biologic behaviors. Inspired by the structure of biological membranes and by internal bio-chemical reactions, Gheorghe Păun initiated the area of membrane computing in 1998 as a theoretical computer science model borrowing many concepts from the cell biology [START_REF] Păun | Computing with Membranes[END_REF]. The model, also called P systems, features a nested membrane-like structure delimiting regions that host objects (modeling cell chemicals) which are transformed locally or communicated to other regions by different types of rules that mimic cell chemical transformations. The evolution of the state of the system is performed by transition steps doing a synchronous parallel execution of all rules. Since a typical model contains hundreds and even thousands of such parallel executions, P systems feature an inherent massive parallelism and the global behavior of the system is emerging from many simple interactions provided by rules application. There exist numerous variants of the model depending on the manipulated objects, used types of rules and the parallel execution strategy.

Membrane Computing is a computing paradigm inspired from the structure and functioning of living cells, and the organization of cells in tissues and other structures, including the brain. It provides distributed parallel computing devices called, generically, P systems. Membrane computing models are an extension of DNA computing, a possible source of novel/useful/interesting computing models. If one wants to model (1) discrete, (2) distributed, (3) parallel, (4) cell-like or tissue-like system, (5) dealing with multisets, [START_REF] Alonso | A Circuit Implementing Massive Parallelism in Transition P Systems[END_REF] evolving by rewriting-like rules, then P systems are obligatory/unavoidable and the unique models dealing with all these features [START_REF] Păun | Membrane Computing: An Introduction[END_REF]. As unconventional computing devices within natural computing, P systems have proved to overcome the well-known limitations imposed by the conventional techniques based on electronic technology. More specifically, the relevance of these computing models can be non-exhaustively summarized in the following points. From a theoretical view, a novel methodology to tackle the famous P versus NP problem has been given [START_REF] Orellana-Martín | P systems with symport/antiport rules: When do the surroundings matter?[END_REF]. The main focus in the area was the introduction of (biologically-inspired) variants of P systems and the further study of their computational power, in particular of their computational completeness, see [START_REF]The Oxford Handbook of Membrane Computing[END_REF] resuming hundreds of articles on this topic. From a practical view, the theoretical investigation led to a series of successful applications in different areas ranging from image processing [START_REF] Díaz-Pernil | Membrane computing and image processing: a short survey[END_REF][START_REF] Yuan | A Resolution-Free Parallel Algorithm for Image Edge Detection within the Framework of Enzymatic Numerical P Systems[END_REF], meta-heuristic algorithms for optimization problems [START_REF] Zhang | A hybrid approach based on differential evolution and tissue membrane systems for solving constrained manufacturing parameter optimization problems[END_REF][START_REF] Zhang | Evolutionary membrane computing: A comprehensive survey and new results[END_REF][START_REF] Zhang | An optimization spiking neural P system for approximately solving combinatorial optimization problems[END_REF], robot controllers [START_REF] Buiu | Development of membrane controllers for mobile robots[END_REF][START_REF] Wang | Design and implementation of membrane controllers for trajectory tracking of nonholonomic wheeled mobile robots[END_REF] and path planning [START_REF] Pérez-Hurtado | A membrane parallel rapidly-exploring random tree algorithm for robotic motion planning[END_REF][START_REF] Pérez-Hurtado | Simulation of rapidly-exploring random trees in membrane computing with P-lingua and automatic programming[END_REF][START_REF] Wang | A modified membrane-inspired algorithm based on particle swarm optimization for mobile robot path planning[END_REF] and fault diagnosis for electrical power systems [START_REF] Rong | Automatic Implementation of Fuzzy Reasoning Spiking Neural P Systems for Diagnosing Faults in Complex Power Systems[END_REF][START_REF] Wang | Fuzzy membrane computing: theory and applications[END_REF][START_REF] Wang | Fault Diagnosis of Electric Power Systems Based on Fuzzy Reasoning Spiking Neural P Systems[END_REF]. Due to its historical background, membrane computing was also used as a modeling framework for biological and ecological subjects such as artificial life [START_REF] Suzuki | Artificial Life Applications of a Class of P Systems: Abstract Rewriting Systems on Multisets[END_REF], photosynthesis [START_REF] Taishin | A Membrane Computing Model of Photosynthesis[END_REF], p53 protein signaling pathway [START_REF] Suzuki | Modeling p53 Signaling Pathways by Using Multiset Processing[END_REF], myxobacterial colony [START_REF] Nash | A P system model of swarming and aggregation in a Myxobacterial colony[END_REF], biopolymer duplication [START_REF] Manca | From biopolymer duplication to membrane duplication and beyond[END_REF]. Please see a recent overview in [START_REF] Zhang | Real-life Applications with Membrane Computing[END_REF].

For twenty years many applications of P systems have arisen, not only in Theoretical Computer Science, but also in Computational Modelling, Robotics, Optimization, etc., as mentioned above. The interest in this area is increasing, given that this unconventional, massively-parallel approach has been demonstrated to provide a powerful, flexible and expressive framework. One example is the 2016 report of the National Research Council of Canada where Membrane Computing appears to be mentioned several times as prominent parts of bio-computing [START_REF] Wiseman | Next Generation Computing[END_REF]. Therefore, there is a need for simulators in order to build model validation tools, assistants for formal verification, environments for virtual experimentation, model calibration tools, etc. So far, most of the simulation software aim at reproducing only the computation of the devices. However, there have been also a research line concerning the implementation of P system parallelism in real parallel platforms. The main hardware employed for this is FPGAs and GPUs, given their efficiency, fast shared memory system and scalability. The objective is twofold: on the one hand, provide efficient simulation tools where the massive, natural P system parallelism do not get serialized and is instead harnessed in order to speed up the simulations; and on the other hand, to explore how this special parallelism can be mapped in current modern computing architectures. It is known that the future of computer architecture will go through non-Von Neumann architectures [START_REF] Basu | iFLEX: A Fully Open-Source High-Density Field-Programmable Gate Array (FPGA)-Based Hardware Co-Processor for Vector Similarity Searching[END_REF]. Natural Computing models can provide solutions where the instructions are closed to, or along with, the data. P systems are just an alternative, like artificial neural networks, for this. We want to shed a light into this alternative from the perspective of efficient implementation, analysing all the challenges and current solutions. In addition, a number of these applications mentioned above only theoretically profit from the potential speed-up promised by the model as there are no truly parallel implementations available. On the other hand, the inherent large scale parallelism of the model has the profound potential for the progress of extreme data processing, especially taking into account that there are not so many widely investigated massively parallel computational paradigms. Thus, an interesting topic is the implementation of membrane computing models on contemporary silicon integrated circuits. This allows to exploit the desirable parallel computational capability of P systems to explore a new orientation for high performance computing (HPC) [START_REF] Katagiri | High-Performance Computing Basics[END_REF][START_REF] Jianping | High performance heterogeneous computing for collaborative visual analysis[END_REF].

As pointed out in [START_REF] Grama | Introduction to Parallel Computing[END_REF] "some unmistakable trends in hardware design indicate that uniprocessor (or implicitly parallel) architectures may not be able to sustain the rate of realizable performance increments in the future". The augmentation of electronic ingredients' density had been subject to the well-known Moore's law for decades. After extraordinary exponential growth of many years, the number of transistor in chips cannot follow this law, at least it cannot be doubled within two years [START_REF]Moore's law organization website[END_REF]. With transistors shrunk to nanoscale, quantum effects stand out [START_REF] Mathur | Beyond the silicon roadmap[END_REF] and the behavior of circuits is not up to expectations. Moreover, even with the further increase of the density, the computational capability growth is not linearly proportional to it [START_REF]The CUDA Handbook: A Comprehensive Guide to GPU Programming[END_REF]. Another knotty problem is the heat dissipation, which would melt the silicon substrate with density level increasing. Traditional semiconductor scaling is predicted to reach an end by about 2024 on the foundation of prior arts [START_REF]IRDS points to chip stacks and new architectures[END_REF]. Parallel computing has the potential to further uplift computing power provided that the density of transistor is constant [START_REF] Trobec | Parallel Computing[END_REF] with multicore and multithread architecture, although heat dissipation and interconnect issues would be challenges [START_REF] Kushwah | Multicore Processor Technology-Advantages and Chanllenges[END_REF].

All these observations give strong arguments in favor of investigation of parallel computing platforms. Before claiming that the era of parallel computing has dawned, one essential question should be clarified: what does parallel computing mean? Though not rigorous, parallel computing implies computing based on the decomposition of the task into a set of concurrently executable operations and the assignment of these operations to multiple parallel processing nodes. The evolution of computer processor scheme, from single-core single central processing unit (CPU) to multiple-core single CPU and multiple-core multiple-CPU frame is an instance of parallel computing. The inherent parallelism of P systems place them in the class of multiple processing nodes computing devices. The mapping of constituents to processing nodes gives rise to different implementation strategies. A common practice in the area of natural computing is to observe natural processes and somehow mimic the corresponding behavior. The way how living cells allocate, organize and coordinate processing nodes has evolved for billions of years. Investigating this magnificent course might help us to handle the multiple cores computing, which has cut a striking figure in the contemporary parallel computing realm.

The large scale distributed parallel processes occurring in vesicle compartments and the vesicle division functionality enlightened from mitosis of living cells are two of the most outstanding advantages of membrane computing. They allow to underlie the foundation for the construction of highly parallel computation platform whose performance, flexibility and scalability outperforms traditional sequential counterparts substantially [START_REF] Zhang | Real-life Applications with Membrane Computing[END_REF]. As a parallel computing paradigm inspired by the structural and functional features of biological membranes, only parallel computing platforms are suitable for the implementation of P systems. More precisely, the limited parallelism of general computers realized by the communication mechanism among the multiple cores of CPU and GPU cannot make full use of the large scale parallelism, non-determinism and other particular attributes that impart an enormous computing potential like creation and dissolution of inner membranes, the self-replication or autopoiesis [START_REF] Uribe | Autopoiesis: The Organization of Living Systems, its Characterization and a Model[END_REF] of the whole cell-like entirety that works as a computing unit, the communication by symport and antiport of objects [START_REF] Alberts | Molecular Biology of the Cell[END_REF], etc. We would like to remark that programming membrane computing algorithms with high level general purpose languages and executing them on the computer represents just a simulation and not a real implementation of P systems [START_REF] Păun | Membrane Computing: An Introduction[END_REF].

Several software-based and hardware-based parallel computing platforms have been developed to implement P systems. The first software-based parallel computing platform was constructed using a cluster of computers [START_REF] Ciobanu | P Systems Running on a Cluster of Computers[END_REF]. This platform achieves good performance and flexibility. Nonetheless, when the size of target P system is increased, the consumption of CPU time and resources caused by the communication between different computers rises dramatically. Moreover, the underlying hardware (a cluster of computers) of this platform cannot be miniaturized closing the way to corresponding membrane computing algorithms to be used in embedded chips and compact controllers which can be employed in robots, automobiles, machine tools, etc. This disadvantage limits the range of applications of such platforms for membrane computing.

Hence, it is important to propose hardware implementations of P systems as specific architectures that do not have the drawbacks related to the traditional ways of implementation. There are two main directions for such research using (1) Field-programmable gate arrays (FPGA) and ( 2) graphical processing units (GPUs) relying on Compute unified device architecture (CUDA) platform. In the first case a completely new parallel circuit is specially designed to implement some variants of P systems. In the second case the pre-defined CUDA parallel platform is used to simulate P systems. The achieved performance and correspondence is smaller in this case, but the development effort is much lower, so finally it becomes an interesting compromise between traditional computers implementations and a highly parallel one using specialized circuits. Also, when implementing membrane computing models in hardware, the main difficulty comes from the fact that there exists a big number of variations of the basic model of P systems having quite distinct characteristics [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF]. This poses a great challenge for the conception of a general computational architecture to implement these various models.

The computation in a P system is a sequence of transitions between configurations. The core problem for implementations is the object distribution problem (ODP) that computes one of the multisets of applicable rules to the current configuration and whose application permits to reach the next configuration. This problem is a particular variant of a more general problem that computes the whole applicable set of multisets of rules for a given configuration and it is known to be NPcomplete [START_REF] Ciobanu | Complexity of evolution in maximum cooperative P systems[END_REF]. Known algorithms and heuristics do not parallelize well, so special heuristics were developed in order to quickly compute the desired multiset of rules. We decided to present these heuristics uniformly in terms of multi-criteria optimization, see Sections 4, 5. Another problem for implementations is that in the general case the model is non-deterministic, so an equitable choice among different possibilities should be provided. However, this is very difficult to achieve and in most of the cases this property is not satisfied.

This paper is organized as follows: Section 2 gives a brief description of the capabilities and of the architecture of hardware used, Section 3 recalls the definition of the most general model of P systems, based on the formal framework [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF]. Next Section 4 expresses object distribution problem (ODP) in terms of multi-criteria optimization and integer linear programming. Section 5 presents an overview of different simulation approaches, including the Direct Non-deterministic Distribution algorithm (DND) algorithm which is the base for handling non-determinism. Next Sections 6, 7 and 8 give more details on existing simulation approaches. At the end, conclusions and future research directions are discussed.

THE SELECTION OF HARDWARE

FPGA hardware

The FPGA is a reconfigurable hardware allowing to prototype digital circuits. The modification of circuits for FPGA is performed by altering the interconnections between circuit elements using a hardware description language, the most common ones being VHSIC Hardware Description Language (VHDL) and Verilog. The design can be performed on several levels of abstraction, ranging from the switch/transistor level until the behavioral level (corresponding to a Mealy machine [START_REF] Mealy | A method for synthesizing sequential circuits[END_REF]).

From structural point of view, an FPGA is an array of configurable logic blocks (CLB) inlaid in the matrix of interconnects. The CLBs work as containers for slice-organized logic cells, which are composed of look-up tables (LUTs) and flip-flops (FFs). The LUTs are used to implement different combinatorial circuits such as basic gates, decoders, encoders and multiplexers. The FFs are the sequential logic components acting as a memory element in the circuit. In order to produce a circuit, the interconnects of CLBs should be reconfigured. The interconnection in FPGA is performed by switch boxes routing signals between its logic blocks. Modern FPGAs use one of the following interconnect technologies: static RAM, flash memory and anti-fuse. The first one dominates the current FPGAs implementations. With the re-programmability, FPGAs comply with the model-oriented hardware implementation quite well.

CUDA-enabled GPU hardware

The multiple computing cores architecture originates from mainframe computers where it was used to improve the clock speed. Nowadays, it is hard to find a PC equipped with only one core. On the other side, the component integration scale of some High-end GPUs has outpaced the CPUs for the booming demand of graphics processing (advanced rendering and 3D vision) [START_REF] Kirk | Programming Massively Parallel Processors: A Hands On Approach[END_REF]. Currently, the GPU is a computing element as powerful as the CPU. Different from FPGAs, there are manufactured parallel architectures in GPUs. The advantage is that developers should just concern about the efficient utilization of these architectures and the drawback is that these frameworks are un-reconfigurable.

Nevertheless, the GPU is not a general processing unit which can handle other computing assignments except for graphics processing. The predicament has changed for the arise of compute unified device architecture, known as CUDA, from the leading chip vendor-NVIDIA corporation. CUDA is a technology that enables general-purpose computing on graphics processing units (GPGPU). Usually, when referring to CUDA, what is referred is not the parallel computing framework but a GPU supporting CUDA. A CUDA-enabled graphics processing unit is an universal parallel computing device which is suitable for the implementing of parallel algorithm models. The parallel computing behavior of CUDA is based on the execution of multiple compute kernels on the GPU. These compute kernels are without physical construction, but based on an abstract parallel programming model. In other words, CUDA does not alter the physical structure of GPU. CUDA programming model is based on heterogeneous computing, where the CPU (host) is the master node that controls the execution flow and launches kernels on the GPU (device) when massive parallelism is required [START_REF] Kirk | Programming Massively Parallel Processors: A Hands On Approach[END_REF]. A kernel is executed by a grid of (thousands of) threads. The grid is a two-level hierarchy, where threads are arranged into thread blocks of equal size. Each block and each thread is unequivocally identified by an identifier. In this way, threads and blocks can be distributed easily to different portions of data, or to compute different instructions. Threads from the same block can be synchronized using barriers, while those belonging to different blocks can only be synchronized by the end of the execution of the kernel.

A GPU contains a global memory, which has the biggest size, but has the longest access time and a shared memory, which is smaller but faster [START_REF] Kirk | Programming Massively Parallel Processors: A Hands On Approach[END_REF]. Although current GPUs contain cache memories, in order to accelerate memory accesses, best performance is achieved when doing it manually. Global memory is accessed by all threads launched in all grids, and also by the host, but shared memory is only accessible by threads in a block. Threads also have fast access to their own registers for single variables, and local memory (which is normally outsourced to global memory). Accesses to memory have to be carefully programmed, so that contiguous portion of data is read by consecutive threads (providing so called coalesced access), since this increases the memory bandwidth utilization.

Nowadays the architecture of GPUs is upgraded to Streaming Multiprocessors (SMs) which are composed of an array of Streaming Processors (SPs), working as computing cores. A thread set consisting of 32 threads named warp is the basic unit which a SM fulfills in its executions. A SM can manage multiple warps which are based on Single-Instruction Multiple-Thread (SIMT) model in effect. Each thread in a warp should commence its processing at the identical program address concurrently, although after beginning, threads can execute independently abiding by a sequential manner. The parallelism of CUDA is terminated when a warp branches or the memory stalls [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF].

Other hardware

There were attempts to simulate certain types of P systems on micro-processor based architectures [START_REF] Gutiérrez | Design of a Hardware Architecture Based on Microcontrollers for the Implementation of Membrane Systems[END_REF]. Although the performance of micro-processors is low, they are economical alternatives suitable for the developing of prototypes and verifying the design methods. The drawback is that in most of the cases the parallel nature of P systems is not exploited at all, often leading to inefficient implementations. It could be interesting to use out-of-order (OoO) execution [START_REF] Smith | Decoupled Access/Execute Computer Architectures[END_REF] processors to tackle this problem, but this possibility was not investigated yet.

Custom application specific integrated circuits (ASIC) can also be employed to implement P systems. However, no such attempts exit at the moment of the writing of this paper. One of the main reasons is that the flexibility of the hardware platform constructed on ASIC is quite insufficient to adapt to the different variants of P systems. However, such attempts could still be meaningful as they would allow to simulate features of P systems on some tailored circuits with different properties, like low power consumption. Another possible direction is to use ASIC for particular commonly occurring computational cores and combine them with software execution [START_REF] Esmaeilzadeh | Neural Acceleration for General-Purpose Approximate Programs[END_REF][START_REF] Venkatesh | QSCORES: Trading dark silicon for scalable energy efficiency with quasi-specific cores[END_REF].

THE MODEL OF P SYSTEMS

We suppose the reader to have a knowledge of basic notions from formal language theory and membrane computing. We refer to [START_REF] Gheorghe | Membrane Computing, Power and Complexity[END_REF][START_REF]The Oxford Handbook of Membrane Computing[END_REF] and [START_REF]Handbook of Formal Languages[END_REF] for missing details. We will also closely follow [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF][START_REF] Verlan | Using the Formal Framework for P Systems[END_REF] for the definitions.

We recall that a multiset can be seen as a set whose elements can have greater than one multiplicity. We will use the string notation for multisets, i.e., a multiset 𝑀 will be represented by a string where the number of occurrences of each letter, corresponds to its multiplicity in 𝑀. We will denote by |𝑀 | the size of the multiset 𝑀 and by |𝑀 | 𝑎 the number of elements 𝑎 in 𝑀.

There exist many variants of P systems, see e.g. [START_REF] Păun | Membrane Computing: An Introduction[END_REF]. In this paper the presentation will be given in terms of the formal framework for P systems introduced in [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF], see also [START_REF] Verlan | Study of Language-Theoretic Computational Paradigms Inspired by Biology[END_REF][START_REF] Verlan | Using the Formal Framework for P Systems[END_REF]. This framework is based on the model of network of cells specifying the structure and rules to be executed as well as on a set of 4 functions giving the semantics of the system. Different combinations of these parameters allow in most of the cases to construct for a given P system a network of cells that strongly bisimulates it (i.e. one step in the one system is simulated by one step in the other one). Moreover, in many cases there is a one-to-one correspondence between rules applied in each system, meaning that they are mostly indistinguishable. Hence, the framework for P systems can be seen as a common language to compare different P systems as well as to express notions related to this area. Moreover, other multiset rewriting-based models (like Petri nets) can be easily expressed into this framework giving a possibility to compare corresponding models. Finally, we would like to remark that in what follows we will use the variant of the framework supposing that the system structure does not change in time. An extension of the framework that permits to take into account notions related to P systems with dynamically evolving structure is given in [START_REF] Freund | A Formalization of Membrane Systems with Dynamically Evolving Structures[END_REF], but corresponding definitions are too complex for the purpose of this paper.

Network of Cells

As pointed out in [START_REF] Freund | Flattening in (Tissue) P Systems[END_REF][START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF][START_REF] Verlan | Study of Language-Theoretic Computational Paradigms Inspired by Biology[END_REF][START_REF] Verlan | Using the Formal Framework for P Systems[END_REF], most types of (static structure) P systems can be seen as variants of parallel multiset rewriting (by using the algorithm called flattening). Since multiset rewriting level is not practical for system description and understanding, a higher-level concept called network of cells was introduced in [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF]. This model augments multiset rewriting with the notion of spatial locations (cells) as well as the corresponding operations and it can be seen as a particular interpretation of the symbols. The works [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF][START_REF] Verlan | Study of Language-Theoretic Computational Paradigms Inspired by Biology[END_REF][START_REF] Verlan | Using the Formal Framework for P Systems[END_REF] define some basic building blocks in terms of network of cells and give several examples of the construction of widespread notions and types of rules in membrane computing using these blocks.

Below, we provide the definition of network of cells, taken from in [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF]. We remark that the definition from [START_REF] Freund | A Formalization of Membrane Systems with Dynamically Evolving Structures[END_REF] is slightly different, however both models coincide when the structure of the system does not evolve.

Definition 3.1 ([35]

). A network of cells of degree 𝑛 ≥ 1 is a construct Π = (𝑛, 𝑉 , 𝑤, 𝐼𝑛𝑓 , 𝑅) where

(1) 𝑛 is the number of cells;

(2) 𝑉 is an alphabet;

(3) 𝑤 = (𝑤 1 , . . . , 𝑤 𝑛 ) where 𝑤 𝑖 ∈ 𝑉 • , for all 1 ≤ 𝑖 ≤ 𝑛, is the finite multiset initially associated to cell 𝑖; (4) 𝐼𝑛𝑓 = (𝐼𝑛𝑓 1 , . . . , 𝐼𝑛𝑓 𝑛 ), where 𝐼𝑛𝑓 𝑖 ⊆ 𝑉 , for all 1 ≤ 𝑖 ≤ 𝑛, is the set of symbols occurring infinitely often in cell 𝑖 (in most of the cases, only one cell, called the environment, will contain symbols occurring with infinite multiplicity); (5) 𝑅 is a finite set of rules of the form

(𝑋 → 𝑌 ; 𝑃, 𝑄)
where 𝑋 = (𝑥 1 , . . . , 𝑥 𝑛 ), 𝑌 = (𝑦 1 , . . . , 𝑦 𝑛 ), 𝑥 𝑖 , 𝑦 𝑖 ∈ 𝑉 • , 1 ≤ 𝑖 ≤ 𝑛, are vectors of multisets over 𝑉 and 𝑃 = (𝑝 1 , . . . , 𝑝 𝑛 ), 𝑄 = (𝑞 1 , . . . , 𝑞 𝑛 ), 𝑝 𝑖 , 𝑞 𝑖 , 1 ≤ 𝑖 ≤ 𝑛 are finite sets of multisets over 𝑉 . We will also use the notation (omitting 𝑝 𝑖 , 𝑞 𝑖 , 𝑥 𝑖 or 𝑦 𝑖 if they are empty)

(1, 𝑥 1 ) . . . (𝑛, 𝑥 𝑛 ) → (1, 𝑦 1 ) . . . (𝑛, 𝑦 𝑛 ) ; [(1, 𝑝 1 ) . . . (1, 𝑝 𝑛 )]; [(1, 𝑞 1 ) . . . (𝑛, 𝑞 𝑛 )].
The above rule is applied as follows: objects 𝑥 𝑖 from cells 𝑖 are rewritten into objects 𝑦 𝑗 produced in cells 𝑗, 1 ≤ 𝑖, 𝑗 ≤ 𝑛, if every cell 𝑘, 1 ≤ 𝑘 ≤ 𝑛, contains all multisets from 𝑝 𝑘 and does not contain any multiset from 𝑞 𝑘 .

By taking for each rule the set of cells that are involved a hypergraph relation, called the structure of the system, is induced. Commonly, tree-like (for P systems) or graph-like (for tissue P systems) relations are considered.

The configuration 𝐶 of Π is defined as an 𝑛-tuple of multisets over 𝑉 (𝑢 1 , . . . , 𝑢 𝑛 ) satisfying 𝑢 𝑖 ∩ 𝐼𝑛𝑓 𝑖 = ∅, 1 ≤ 𝑖 ≤ 𝑛.

In order to define the computation in network of cells according to some derivation mode 𝛿 the following functions should be specified:

• 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, C, 𝛿) -the function taking a system Π, a configuration C and a derivation mode 𝛿 and yielding the set of multisets of rules of Π that can be applied to C. • 𝐴𝑝𝑝𝑙𝑦 (Π, C, 𝑅) -the function allowing to compute the configuration obtained by the parallel application of the multiset of rules 𝑅 to the configuration C. • 𝐻𝑎𝑙𝑡 (Π, C, 𝛿) -a predicate that yields true if C is a halting configuration of the system Π (in some derivation mode 𝛿). • 𝑅𝑒𝑠𝑢𝑙𝑡 (Π, C) -a function giving the result of the computation of the P system Π when the halting configuration C has been reached. Then the computation is a sequence of transitions where each transition step 𝐶 ⇒ 𝐶 ′ is defined as 𝐶 ′ = 𝐴𝑝𝑝𝑙𝑦 (Π, 𝐶, 𝑅), for some 𝑅 ∈ 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿). As usual, this sequence starts with the initial configuration and ends with the final configuration for which the halting predicate 𝐻𝑎𝑙𝑡 yields true. In more formal terms, the result of the computation of a network of cells Π = (𝑛, 𝑉 , 𝑤, 𝐼𝑛𝑓 , 𝑅) working in the derivation mode 𝛿 is defined as follows (we refer to [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF] for more technical details):

𝑅𝑒𝑠𝑢𝑙𝑡 (Π) = {𝑅𝑒𝑠𝑢𝑙𝑡 (Π, 𝑧) : 𝑤 ⇒ * 𝑧, 𝐻𝑎𝑙𝑡 (Π, 𝑧, 𝛿) = 𝑡𝑟𝑢𝑒 and 𝐻𝑎𝑙𝑡 (Π, 𝑥, 𝛿) = 𝑓 𝑎𝑙𝑠𝑒 for any 𝑥 : 𝑤 ⇒ * 𝑥 ⇒ + 𝑧}
We remark that [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF] gives an algorithm to compute the set of multisets of applicable rules, denoted as 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛), the algorithm to compute 𝐴𝑝𝑝𝑙𝑦 (Π, C, 𝑅) and several definitions for 𝐻𝑎𝑙𝑡 and 𝑅𝑒𝑠𝑢𝑙𝑡 functions. The most common way of halting is the total halting, which means that there are no more applicable rules and the most common way of getting the result is to consider the multiset of objects present in the halting configuration at some predefined cell.

At each step, the set of multisets of applicable rules 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛) can be restricted using a derivation mode, denoted by 𝛿, which specifies which sub-multisets are chosen for the next step application. The most common example is the maximally parallel derivation mode (𝑚𝑎𝑥) which is defined as follows:

𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) = {𝑅 ⊆ 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛) | 𝑅 ′ ∈ 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛) : 𝑅 ′ ⊋ 𝑅}.
The filter condition states that only non-extensible multisets are considered in this derivation mode. We remark that there might be several such multisets, hence for the application a nondeterministic choice is used to take one of them. We also refer to [START_REF] Verlan | Using the Formal Framework for P Systems[END_REF] for a description of several other derivation modes.

Example 3.2. Consider the system Π = (𝑂, 𝑤 1 , 𝑅), having the alphabet 𝑂 = {𝑎, 𝑏, 𝑐} and the set of rules 𝑅 = {𝑟 1 :

(1, 𝑎𝑏) → (1, 𝑎𝑏𝑐); 𝑟 2 : (1, 𝑏𝑏𝑐) → (1, 𝑎𝑏𝑏)}. Consider the configuration 𝐶 = (1, 𝑎 3 𝑏 4 𝑐 2 ). Then, 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛) = 𝑟 1 , 𝑟 2 1 , 𝑟 3 1 , 𝑟 2 , 𝑟 2 2 , 𝑟 1 𝑟 2 , 𝑟 2 1 𝑟 2 .
The maximally-parallel set of multisets of rules is the following:

𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) = 𝑟 3 1 , 𝑟 2 2 , 𝑟 2 1 𝑟 2 The application of the multiset of rules 𝑟 2 1 𝑟 2 on 𝐶 yields (1, 𝑎 4 𝑏 4 𝑐 3 ). Example 3.3. Consider the system Π = (𝑂, 𝑤 1 , 𝑤 2 , 𝑅), with 𝑂 = {𝑎, 𝑏, 𝑐} and 𝑅 = {𝑟 1 : (1, 𝑎) (2, 𝑏) → (1, 𝑏) (2, 𝑎); 𝑟 2 : (1, 𝑏) → (2, 𝑏)}. Consider the configuration 𝐶 = (1, 𝑎 4 ) (2, 𝑏).
It is easy to observe that at each step only a single rule is applicable (𝑟 1 or 𝑟 2 alternatively). The only possible sequence of rule applications is (𝑟 1 𝑟 2 ) 4 after which no rule is applicable anymore. This sequence moves all symbols 𝑎 from cell 1 to cell 2. We remark that the above rules do not rewrite objects, but only move them in the structure.

Examples of P systems

Here we will give the formal framework description of some main variants of P systems that were targeted for an implementation.

Transitional (cell-like) P systems [START_REF] Orellana-Martín | Minimal cooperation as a way to achieve the efficiency in cell-like membrane systems[END_REF]. This model considers that cells are organized in a tree structure and uses rules of following type: (𝑖, 𝑢) → (𝑖, 𝑢 ′ ) ( 𝑗, 𝑢 ′′ ) (𝑘 1 , 𝑣 1 ) . . . (𝑘 𝑚 , 𝑣 𝑚 ), where 𝑗 is the parent of 𝑖 and 𝑖 is the parent of 𝑘 1 , . . . , 𝑘 𝑚 , 𝑚 ≥ 0. When |𝑢 | > 1 corresponding rules are called cooperative.

Symport/antiport P systems. This variant considers that objects are not transformed but rather moved in a tree-like or graph-like structure. This corresponds to rules of form (𝑖, 𝑢) ( 𝑗, 𝑣) → (𝑖, 𝑣) ( 𝑗, 𝑢) or (𝑖, 𝑢) → ( 𝑗, 𝑢).

Populational Dynamics P (PDP) systems [START_REF] Bernardini | Population P Systems[END_REF][START_REF] Sánchez-Karhunen | Modelling complex market interactions using PDP systems[END_REF]. The structure of this model corresponds to several trees that have their roots linked together. There are two types of rules: (1) working inside some tree: (𝑖, 𝑢) ( 𝑗, 𝑣) → (𝑖, 𝑢 ′ ) ( 𝑗, 𝑣 ′ ), where 𝑖 is parent of 𝑗 and (2) communication between tree roots: (𝑟 𝑖 , 𝑥) → (𝑟 1 , 𝑥 1 ) . . . (𝑟 𝑘 , 𝑥 𝑘 ), where 𝑟 1 , . . . , 𝑟 𝑘 are the numbers of the corresponding tree roots. Moreover, each rule has a probability associated to it, so the derivation mode is maximally parallel, followed by a probabilistic choice between rules having the same left-hand-side.

Spiking Neural P systems [START_REF] Ren | Generating context-free languages using spiking neural P systems with structural plasticity[END_REF][START_REF] Jimenez | Matrix representation and simulation algorithm of spiking neural P systems with structural plasticity[END_REF]. This model has a graph structure and restricts the alphabet of the system to be a single letter, however permitting and forbidding conditions are replaced by a regular expression check. The rules are of form (𝑖, 𝑎 𝑘 ) → (𝑘 1 , 𝑎 𝑥 1 ) . . . (𝑘 𝑚 , 𝑎 𝑥 𝑚 ); (𝑖, 𝐸), where 𝐸 is the regular expression checking for the contents of cell 𝑖, 𝑎 is the single symbol used in the alphabet and 𝑘 1 , . . . , 𝑘 𝑚 are the cells linked to cell 𝑖. The derivation mode is sequential at the level of each cell (only one rule per cell is applied) and maximally parallel at the level of all cells.

(Enzymatic) Numerical P Systems [START_REF] Pavel | Enzymatic numerical P systems -a new class of membrane computing systems[END_REF][START_REF] Păun | Membrane Computing and Economics: Numerical P Systems[END_REF]. A multiset 𝑀 can also be seen as a function 𝑀 : 𝑉 → N. The model of numerical P systems extends the notion of multiset to 𝑀 : 𝑉 → R 𝑛 . Hence, the configuration is a vector of real-valued variables and rules are equations describing how to update each variable during each step. [START_REF] Gazdag | A new method to simulate restricted variants of polarizationless P systems with active membranes[END_REF]. This model corresponds to transitional P systems with an additional rule allowing to create new child cells.

P Systems with Active Membranes

REDUCTION OF THE COMPUTATION STEP TO MULTI-CRITERIA OPTIMIZATION

The computation in P systems corresponds to a finite sequence of transitions between configurations. Each transition 𝐶 ⇒ 𝐶 ′ can be seen as a sequence of 4 steps: (1) computing the set of multisets of applicable rules (𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛)), (2) restricting it according to the derivation mode 𝛿 (𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿)), (3) choosing one element from it 𝑅 and (4) applying 𝑅 to the configuration 𝐶 yielding 𝐶 ′ = 𝐴𝑝𝑝𝑙𝑦 (Π, 𝐶, 𝑅). The first step is the most time consuming, as the corresponding problem is generally NP-hard. As we show below that in some cases it is possible to reduce the first 3 steps (the computation of an element from 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿)) to a multi-criteria optimization problem. This will allow us to express different algorithms used to compute this step by different authors in a common language for a better comparison. Another theoretical advantage of such reduction is an extensive standard vocabulary and a plethora of solving methods existing in the optimization area, that could be possibly applied for large-scale P systems. We would like to remark one more time that the presented reductions are not possible in the general case, however they are possible for many concrete cases.

Preliminaries

A multi-criteria optimization problem (MCOP) is an optimization problem that involves multiple objective functions. In mathematical terms it can be stated as max(𝑓 1 (𝑥), 𝑓 2 (𝑥), . . . , 𝑓 𝑚 (𝑥))

subject to 𝑥 ∈ 𝑋
where 𝑚 ≥ 2 and 𝑋 is the set of feasible vectors (or solutions). This set is usually defined by some constraint functions. We can also consider the objective function as a vector: 𝑓 : 𝑋 → R 𝑚 , 𝑓 = (𝑓 1 , . . . , 𝑓 𝑚 ). For a feasible solution 𝑥, the vector 𝑧 = 𝑓 (𝑥) is called an objective vector or an outcome [START_REF] Kaliszewski | Multiple Criteria Decision Making by Multiobjective Optimization-A Toolbox[END_REF][START_REF] Legriel | Multi-Criteria Optimization and its Application to Multi-Processor Embedded Systems[END_REF].

When corresponding functions as well as 𝑓 𝑖 , 1 ≤ 𝑖 ≤ 𝑚 are linear, we speak about a multicriteria linear optimization problem. We also remark that as for classical optimization problems the objective functions are minimized; the other cases like maximization or hybrid min/max can be easily reduced to the minimization one. When further 𝑋 ⊆ N 𝑘 , 𝑘 > 0 and 𝑓 : 𝑋 → N 𝑚 we speak about an integer multi-criteria linear optimization problem (IMCLOP).

In multi-criteria optimization, typically there is no feasible solution that minimizes all objective functions simultaneously. Hence, the main attention is focused on solutions that cannot be improved in any of the objectives without degrading some other objective(s). Such solutions are called Paretooptimal. Formally, they are defined as preimages of maximal elements of the outcomes, which are also called Pareto front. Definition 4.1. A vector 𝑥 ∈ 𝑋 is Pareto-optimal for a MCOP (defined as above) iff there is no other vector 𝑦 ∈ 𝑋 for which 𝑓 (𝑥) < 𝑓 (𝑦), where 𝑢 < 𝑣 iff 𝑢 𝑖 ≤ 𝑣 𝑖 , 1 ≤ 𝑖 ≤ 𝑘, and ∃𝑗, 1 ≤ 𝑗 ≤ 𝑚 such that 𝑢 𝑗 < 𝑣 𝑗 .

Example 4.2. Consider the following problem:

max(𝑟 1 , 𝑟 2 ) subject to 𝑟 1 ≤ 5 𝑟 1 + 2𝑟 2 ≤ 6 𝑟 2 ≤ 3 𝑟 1 ∈ N, 𝑟 2 ∈ N
The corresponding feasible solutions and Pareto-optimal solutions are shown in Fig. 1. We note that Pareto-optimal solutions are (5,0), (4,1), (2,2) and (0,3).

Rule choice as IMCLOP

It is not difficult to see that the problem of the computation of elements from 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) can be reduced to IMCLOP. For the first time it was noticed in [START_REF] Alhazov | Maximally Parallel Multiset-Rewriting Systems: Browsing the Configurations[END_REF], but without any further development. For simplicity, we consider that Π is a transitional P system and has only one membrane (and no environment). If this is not the case, we can apply the flattening procedure reducing it to one membrane [START_REF] Freund | Flattening in (Tissue) P Systems[END_REF][START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF]. So Π = (𝑂, 𝑤 1 , 𝑅).

Let 𝑅 = {𝑟 1 , . . . , 𝑟 𝑚 } and 𝑂 = {𝑎 1 , . . . , 𝑎 𝑛 }. Consider that 𝑟 𝑖 : (1, 𝑢 𝑖 ) → (1, 𝑣 𝑖 ), 1 ≤ 𝑖 ≤ 𝑚.
Let 𝐶 be the current configuration and let 𝐶 𝑎 = |𝐶 | 𝑎 , 𝑎 ∈ 𝑂.

We will consider a set of variables 𝑥 𝑖 , 1 ≤ 𝑖 ≤ 𝑚, that indicate the cardinality of corresponding rules in some parallel multiset M ∈ 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑎𝑠𝑦𝑛), |M | 𝑟 𝑖 = 𝑥 𝑖 , 1 ≤ 𝑖 ≤ 𝑚. Then the feasible set 𝑋 of (asynchronous) solutions is defined by the following inequalities:

𝑚 𝑖=1 |𝑢 𝑖 | 𝑎 𝑥 𝑖 ≤ 𝐶 𝑎 , ∀𝑎 ∈ 𝑂, (1a)
𝑥 𝑖 ∈ N, 1 ≤ 𝑖 ≤ 𝑚, (1b)
Inequalities (1a) state that the sum of all consumed objects is included in 𝐶. Technically, for each object 𝑎 ∈ 𝑂 it is verified that the weighted sum (by the number of symbols 𝑎 in the left-hand-side) of rule cardinalities is smaller or equal than the number of objects 𝑎 in 𝐶. We remark that system (1) can be also seen as a system of Diophantine equations and corresponding solutions are exactly describing the set 𝑋 of feasible solutions.

The IMCLOP corresponding to the computation of 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) is defined by:

max(𝑥 1 , . . . , 𝑥 𝑚 ) (2)
subject to (𝑥 1 , . . . , 𝑥 𝑚 ) ∈ 𝑋 It should be clear that Pareto-optimal solutions represent exactly the multiplicities of rules for some maximally parallel solution. So this system has 4 Pareto-optimal solutions: (5, 0), (4, 1), (3, 3) and (0, 3), corresponding to multisets of rules 𝑟 5 1 , 𝑟 4 1 𝑟 2 , 𝑟 2 1 𝑟 2 2 and 𝑟 3 2 , which are exactly the maximal multisets of rules applicable to 𝐶.

Tentative solutions

One of the traditional approaches to solve multi-criteria optimization problems is called scalarization. It consists in reducing the corresponding MCOP to a single objective optimization problem by using a real-valued scalarizing function involving the objective functions and additional scalar or vector parameters and variables. This can also imply additional restrictions to the feasible set based on the newly introduced variables.

One of the "simplest" methods to solve multi-criteria problems is the weighted sum method, where we solve the following single objective optimization problem max

𝑥 ∈𝑋 𝑚 𝑘=1 𝜆 𝑘 𝑓 𝑘 (𝑥) (3)
The weighted sum problem (3) is constructed using the scalar vector product of the vector of objective functions 𝑓 and the vector of non-negative weights 𝜆 ∈ R 𝑚 as a parameter. It is known that it allows to compute all Pareto-optimal solutions for convex problems by varying 𝜆 [START_REF] Ehrgott | Multicriteria Optimization[END_REF].

In the literature on P systems some simple variants of the weighted sum method can be found. In [START_REF] Agrigoroaiei | Evolving by Maximizing the Number of Rules: Complexity Study[END_REF][START_REF] Reina-Molina | Integer Linear Programming for Tissuelike P Systems[END_REF] the vector 𝜆 = (1, . . . , 1) is considered (so the objective function is the sum of all variables). However, in this case only maximally parallel rulesets having a maximal number of rules are obtained. In terms of the formal framework [START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF] this corresponds to 𝑚𝑎𝑥 𝑟𝑢𝑙𝑒𝑠 𝑚𝑎𝑥 derivation mode.

Another attempt was done in [START_REF] Ciobanu | New Strategies of Using the Rules of a P System in a Maximal Way: Power and Complexity[END_REF][START_REF] Ciobanu | Complexity of evolution in maximum cooperative P systems[END_REF] where the parameters 𝜆 𝑘 correspond to the size of the left-hand-side of rules (𝜆 𝑘 = |𝑢 𝑘 |, 𝑟 𝑘 : 𝑢 𝑘 → 𝑣 𝑘 ). Such optimization problem finds only maximally parallel solutions involving the maximal number of objects, corresponding to 𝑚𝑎𝑥 𝑜𝑏 𝑗𝑒𝑐𝑡𝑠 𝑚𝑎𝑥 mode in terms of the formal framework.

In [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF] the set of maximally parallel multisets of rules is expressed as solutions of a system of Diophantine equations (roughly equations (1a)) with an additional constraint to be satisfied on a solution, expressed as another system of Diophantine equations.

Finally, in [START_REF] Arteta | Algorithm for Application of Evolution Rules Based on Linear Diofantic Equations[END_REF] the set of maximally parallel multisets of rules can be expressed as solutions to a system of equations defining some Diophantine sets. While the construction is similar to the one we give below, it is not trivial to manipulate Diophantine sets and it is not clear how to express the constraints as a single system of equations. Below, we show how to handle this problem by using same construction as for handling multiple "either-or"constraints in integer linear programming.

As above we consider that Π has only one membrane (and no environment). So, Π = (𝑂, 𝑤 1 , 𝑅).

Let 𝑅 = {𝑟 1 , . . . , 𝑟 𝑚 } and 𝑂 = {𝑎 1 , . . . , 𝑎 𝑛 }. Consider that 𝑟 𝑖 : (1, 𝑢 𝑖 ) → (1, 𝑣 𝑖 ), 1 ≤ 𝑖 ≤ 𝑚. Let 𝐶 be the current configuration and let 𝐶 𝑎 = |𝐶 | 𝑎 , 𝑎 ∈ 𝑂.
In addition to the 𝜆 vector above we introduce an integer parameter 𝑀 ∈ N having a value that is sufficiently big (in fact it should be greater than the maximal possible value of any variable 𝑥 𝑖 , 1 ≤ 𝑖 ≤ 𝑚 of any feasible solution). We remark that in general case of maximal parallelism it is impossible to limit the maximal values of variables, however, for practical applications, it is often possible to find such a bound.

We will modify the system (1) to construct a system of inequalities whose integer solutions will be maximally parallel multisets of rules for the configuration 𝐶.

𝑚 𝑖=1 |𝑢 𝑖 | 𝑎 𝑥 𝑖 ≤ 𝐶 𝑎 , ∀𝑎 ∈ 𝑂, ( 4a
) 𝑚 𝑖=1 |𝑢 𝑖 | 𝑎 𝑥 𝑖 + |𝑢 𝑘 | 𝑎 + 𝑀𝑧 𝑘 𝑎 ≥ 𝐶 𝑎 + 1, 1 ≤ 𝑘 ≤ 𝑚, 𝑎 ∈ 𝑂, |𝑢 𝑘 | 𝑎 > 0, ( 4b
) 𝑎 ∈𝑂 𝑧 𝑖 𝑎 = 𝑁 𝑖 -1, 1 ≤ 𝑖 ≤ 𝑚, 𝑁 𝑖 = 𝑎 ∈𝑂 sgn(𝑧 𝑖 𝑎 ) (4c) 𝑥 𝑖 ∈ N, 1 ≤ 𝑖 ≤ 𝑚, (4d)
𝑧 𝑖 𝑎 ∈ {0, 1}, 1 ≤ 𝑖 ≤ 𝑚, 𝑎 ∈ 𝑂. (4e)
Inequalities (4a) are the same as (1a) and they state that the sum of all consumed objects is included in 𝐶.

Inequalities (4b) state the maximality property of the rule set defined by 𝑥 1 , . . . , 𝑥 𝑚 . It verifies that for each rule there exist at least one object whose remaining quantity is not sufficient to apply this rule. They are based on multiple "either-or" constraints representation in ILP. The big value of 𝑀 and the inequalities (4c) ensure that only one constraint from (4b) will be considered (the other ones will be satisfied because of the big value of 𝑀).

Hence, any solution 𝑥 1 , . . . , 𝑥 𝑚 satisfying the system of inequalities (4) corresponds to a Paretooptimal solution of (3), hence to a maximally parallel rule set M = 𝑟 𝑥 1 1 . . . 𝑟 𝑥 𝑚 𝑚 applicable to configuration 𝐶. We also remark that from the construction given above, it immediately follows that system (4) is Diophantine. ). Then, we construct an ILP according to the rules above (we recall that 𝑀 is a big integer number): Derived from (4a):

𝑥 1 + 𝑥 2 ≤ 2 𝑥 1 + 𝑥 3 ≤ 3 𝑥 1 ≤ 2
Derived from (4b) for 𝑟 1 :

𝑥 1 + 𝑥 2 + 1 + 𝑀𝑧 1 𝑎 ≥ 3 𝑥 1 + 𝑥 3 + 1 + 𝑀𝑧 1 𝑏 ≥ 4 𝑥 1 + 1 + 𝑀𝑧 1 𝑐 ≥ 3 Derived from (4c) for 𝑟 1 : 𝑧 1 𝑎 + 𝑧 1 𝑏 + 𝑧 1 𝑐 = 2 Derived from (4b) for 𝑟 3 : 𝑥 1 + 𝑥 2 + 1 ≥ 3
Derived from (4b) for 𝑟 3 :

𝑥 1 + 𝑥 3 + 1 ≥ 4
It is not difficult to see that this system has 3 solutions: (2, 0, 1), (1, 1, 2) and (0, 2, 3), corresponding to multisets of rules 𝑟 2 1 𝑟 3 , 𝑟 1 𝑟 2 𝑟 2 3 and 𝑟 2 2 𝑟 3 3 , which are exactly the maximal multisets of rules applicable to 𝐶.

We remark that all solutions of (4),(3) might be tedious to obtain. For the simulation purposes, only one such solution is necessary.

As previously mentioned, using the weighted sum scalarization technique it is theoretically possible to reach any single point from the Pareto front by choosing appropriate values of the parameter vector 𝜆. However, by using equations (4b) it becomes much easier to chose corresponding parameters as the feasible set is restricted only to Pareto-optimal values. For example, using max 𝑥 1 as objective function with the constraints from Example 4.4, would yield the solution (2, 0, 1). By using min(𝑥 1 + 𝑥 2 + 𝑥 3 -5) as the objective function, the solution (0, 2, 3) is obtained (the above constraint allows to consider only multisets of rules of size 5).

In the literature on P systems, there are other examples of the construction of a single Paretooptimal solution having certain properties. One of them is the DND algorithm introduced in [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF] for the FPGA simulations that we discuss below.

In what follows, we explain briefly the functioning of the DND algorithm in the view of equations [START_REF]IRDS points to chip stacks and new architectures[END_REF]. First a random rule permutation is computed (corresponding to a random permutation of variable indices). Next, during the forward step a random value (bounded by the number of possible applications) for the number of each rule applications is taken. This corresponds to finding the values of 𝑥 𝑖 , satisfying the constraints (1a). Finally, during the backward step, the frequency of each rule is increased until it cannot be applied anymore. This step corresponds to the elimination of the dominated, i.e. smaller, solutions for the system (1) yielding only Pareto-optimal ones. The described procedure can also be seen as a combination of scalarizing and of the lexicographic method, which is a method from the family of a priori methods for multi-criteria optimization. More precisely, it corresponds to a solution of a series of optimization problems, each of them bounded by the parameter corresponding to the choice of the maximal rule multiplicity for the forward step. The backward step also corresponds to a series of optimization problems that just reach the maximum for each component (like in the lexicographical method).

A simpler variant of the DND algorithm that does not perform the initial rule permutation can be found in [START_REF] Victor | Hardware Implementation of a Bounded Algorithm for Application of Rules in a Transition P-System[END_REF][START_REF] Tejedor | Algorithm of rules applications based on competitiveness of evolution rules[END_REF][START_REF] Tejedor | Optimizing Evolution Rules Application and Communication Times in Membrane Systems Implementation[END_REF]. Other variations of DND algorithm were used for CUDA-based simulations and are discussed later. A different approach is used in [START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF]. It supposes that for a P system Π working in the derivation mode 𝛿 there exists an easily computable function 𝑁 𝐵𝑉 𝑎𝑟𝑖𝑎𝑛𝑡𝑠 (Π, 𝐶, 𝛿) that for any configuration 𝐶 gives the number of solutions of (4). Next, it also supposes that there exists an easily computable function 𝑉 𝑎𝑟𝑖𝑎𝑛𝑡 (Π, 𝐶, 𝛿, 𝑛) that for each integer 𝑛 (up to the corresponding value) yields the corresponding solution (the used method is similar to the decoding of a number in the combinatorial number system). It is clear that such functions do not exist for any variant of P system. The papers [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF][START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF] give some sufficient relations between rules that allow to construct the above functions for a P system working in set-maximal (𝑠𝑚𝑎𝑥) derivation mode.

IMPLEMENTATION OF P SYSTEMS USING HARDWARE

The discussion about the hardware implementation of membrane computing concerns mainly the simulation of concrete variants (sometimes even examples) of P systems using a dedicated hardware (for ASIC and FPGA) or firmware (for CUDA). In most of the cases a (single) particular class of P systems is simulated. The simulator is composed from two parts: (1) the hardware simulator for a concrete system and (2) the software generator of hardware simulators, which based on input parameters (rules, membranes, initial configuration etc.) generates the code for the dedicated hardware simulator. In this section we will mainly discuss the structure of the corresponding hardware simulators, as the generator part is more-or-less following standard compiler construction techniques. We will concentrate on three points, which are the most important for a hardware implementation: (a) the representation of the configuration (multisets of objects and the membranes), (b) the representation of rules and their parallel application and (c) handling of the non-determinism.

Data organization

Membranes. We remark that there are no compartments in silicon circuits, hence the notion of membranes is relatively difficult to represent directly. Nevertheless, according to [START_REF] Păun | Membrane Computing: An Introduction[END_REF], a membrane is just an idealized concept without internal structures. The main functionality of membranes is to perform a topological division of the space allowing P systems to compute in a distributed manner (based on an correspondence between the membrane and its contents). Hence, the spatial placement and size of membrane are not important, only the inter-relationship among them matters. Moreover, it is known that any membrane structure can be reduced (flattened) to just a single membrane, see [START_REF] Freund | Flattening in (Tissue) P Systems[END_REF][START_REF] Freund | A Formal Framework for Static (Tissue) P Systems[END_REF] for more details. The existing hardware simulators adapt in most cases this last point of view, where the membrane structure is not physically implemented on the device. As exceptions from this rule we cite [START_REF] Petreska | A Reconfigurable Hardware Membrane System[END_REF] and [START_REF] Van Nguyen | A Region-Oriented Hardware Implementation for Membrane Computing Applications[END_REF] (region-based), which implicitly implement the membrane topology by using dedicated buses and message passing in the corresponding circuits.

Configuration. The representation of the configuration in all cases is done as a vector of nonnegative integers (stored in the memory/registers of the device). We remark that this vector corresponds to a flattened system, so it is relatively big and sparse, as it encodes each object/membrane pair. For performance reasons, it is physically split in several places to be closer to the processing units (as routing is relatively expensive). In the case of [START_REF] Van Nguyen | A Region-Oriented Hardware Implementation for Membrane Computing Applications[END_REF][START_REF] Petreska | A Reconfigurable Hardware Membrane System[END_REF] it can be argued that corresponding parts are internalized into the corresponding region circuit, as the access to corresponding values is not direct and it is done by message passing. In CUDA, only the portions of the array representing the configuration in the corresponding membranes to be processed are loaded by the cores.

Evolution Rules. As it can be seen from definition 3.1 in the simple case (without permitting and forbidding) rules can be defined by 2 integer vectors indicating the multiplicities of corresponding objects in the left-hand-side and right-hand-side of each rule. This gives a natural rule representation as two vectors stored in the memory/registers. Then a specific circuit/module verifies the applicability of rules and performs their application. Most variants of hardware implementations use this idea, however in [START_REF] Van Nguyen | Balancing Performance, Flexibility, and Scalability in a Parallel Computing Platform for Membrane Computing Applications[END_REF][START_REF] Van Nguyen | An Implementation of Membrane Computing Using Reconfigurable Hardware[END_REF][START_REF] Petreska | A Reconfigurable Hardware Membrane System[END_REF] each rule is encoded in hardware as a specific circuit that verifies the needed resources and performs the rule application. In the case of P systems with active membranes, there might be several cells with the same label, so the corresponding CUDA simulator uses a split representation: 2 integer vectors, as above, for all rules, and an index array with one position per rule giving for each cell the rules it contains.

Object Distribution Problem and Non-Determinism

The hardware implementation of P systems faces the problem of the computation of the applicable rule set according to some derivation mode (usually maximally parallel). More precisely, an efficient way to compute and represent an element from 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿) is required [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF][START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF]. The difficulty of the problem is that rules can compete for same objects, so increasing the number of occurrences for one rule, may decrease the application possibilities for another one. Another important problem is to ensure that a non-deterministic choice among all possibilities is performed. In [START_REF] Arroyo | A Hierarchical Architecture with Parallel Comunication for Implementing P Systems[END_REF][START_REF] Gutiérrez | Hardware and Software Architecture for Implementing Membrane Systems: A Case of Study to Transition P Systems[END_REF][START_REF] Martínez | Hardware Circuit for the Application of Evolution Rules in a Transition P-system[END_REF][START_REF] Victor | Hardware Implementation of a Bounded Algorithm for Application of Rules in a Transition P-System[END_REF], hardware architectures aiming at parallel processing and communication, and the application of rules are developed. In [START_REF] Angel | Towards an Electronic Implementation of Membrane Computing: A Formal Description of Non-deterministic Evolution in Transition P Systems[END_REF], a formal exposition of non-deterministic evolution in transition P systems was suggested.

We call the first problem as object distribution problem (ODP). It consists in the computation of the set 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿) (or of an element from this set). As discussed in Section 4 in terms of multi-criteria optimization this corresponds to the computation of the corresponding Pareto front (or an element of it).

In [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF] different algorithms solving ODP are classified in direct and indirect ones. In the direct approach, the corresponding multiset is directly constructed by the algorithm. In terms of MCOP this corresponds to a particular fixed scalarization. The indirect approaches are based on the observation that the solution number is finite because the solution space is bounded by the size of the configuration. Hence, a heuristic or brute-force approach can be used to explore this bounded space. However, since it is an overestimation, there might be visited elements that are not valid solutions. Hence, the algorithms are iterative and explore the whole space until a valid solution is encountered. In terms of MCOP this corresponds to different searches through the space limited only by the maximal values for each axis.

Sometimes it is not easy to classify an algorithm in one of these categories. We will classify an algorithm as a direct approach if its main goal is to construct a valid multiset of rules. Otherwise, if an algorithm is exploring different solutions until it reaches a valid one, it will be classified as indirect. We will use this classification to overview different strategies for ODP solution known in the literature.

Indirect approaches. Generally, the enumeration of all possible solutions and their verification one by one until a correct solution is obtained is the simplest method for the indirect approach [START_REF] Fernández | Massively Parallel Algorithm for Evolution Rules Application in Transition P Systems[END_REF]. Before the first correct solution is obtained, some invalid solutions should be rejected. This approach is called indirect straightforward approach [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF]. Taking into account that it is not viable to enumerate all possible solutions for many problems, the feasibility of the approach is low. However, the performance of the algorithm suggests its use as to compute the floor values for the object distribution problem. Another indirect approach discussed in [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF][START_REF] Van Nguyen | A Region-Oriented Hardware Implementation for Membrane Computing Applications[END_REF] called indirect incremental approach investigates a strategy generating possible solutions in rounds. Other attempts based on a similar idea but with different rule elimination strategies were done in [START_REF] Fernández | Decision Trees for Applicability of Evolution Rules in Transition P Systems[END_REF][START_REF] Gil | Delimited Massively Parallel Algorithm based on Rules Elimination for Application of Active Rules in Transition P Systems[END_REF][START_REF] Gutiérrez | Optimizing Membrane System Implementation with Multisets and Evolution Rules Compression[END_REF][START_REF] Francisco | Delimited Massively Parallel Algorithm Based on Rules Elimination for Application of Active Rules in Transition P Systems[END_REF][START_REF] Tejedor | Algorithm of rules applications based on competitiveness of evolution rules[END_REF][START_REF] Tejedor | Algorithm of Active Rule Elimination for Application of Evolution Rules[END_REF] Direct approaches. In contrast to indirect approaches, the direct approach fabricates a solution straightforwardly rather than identifying a number of possible solutions before a solution is confirmed.

The simplest approach is the direct straightforward approach. As defined in the paper [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF], in this approach "all the solutions to the object distribution problem are given as input, and one of these solutions is simply selected at random". While in the same paper it is argued that such approach is infeasible for an arbitrary configuration and rule types, it can still be applied in a big number of cases. As shown in [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF][START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF], if at each step the number of solutions can be expressed as a the number of words of some length in a regular language, then it becomes possible to compute the solution only based on its number. In [START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF] it is shown that the corresponding class of P systems is quite large and also that this method is particularly interesting for bounded derivation modes like the set-maximal derivation mode (called also flat mode) where the rules are chosen in a set-maximal way (instead of the multiset maximal way).

Another variant of the direct approach is the Direct Non-deterministic Distribution algorithm (DND) proposed in [START_REF] Van Nguyen | An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware[END_REF]. A similar algorithm can also be found in [START_REF] Gil | Parallel algorithm for P Systems implementation in multiprocessors[END_REF][START_REF] Gil | Fast Linear Algorithm for Active Rules Application in Transition P Systems[END_REF]. This algorithm works in 2 phases. At the first phase all rules (initially randomly shuffled) except one are selected to be applied a random number of times below its maximal applicability value. In the second phase, all the rules are taken in the converse order (with the first one being the rule excluded at the previous step) and their applicability is increased up to the maximal still possible value. A variant of DND, named DND-P, became popular in the simulation of Population Dynamics P (PDP) systems [START_REF] Miguel | A new simulation algorithm for multienvironment probabilistic P systems[END_REF]. Together with another algorithm Direct distribution based on Consistent Blocks Algorithm (DCBA) [START_REF] Miguel | DCBA: Simulating Population Dynamics P Systems with Proportional Object Distribution[END_REF] it was employed for the engine of the PDP system simulator on CUDA [START_REF] Miguel | Parallel simulation of Population Dynamics P systems: updates and roadmap[END_REF].

Non-determinism. One of the difficulties of the above approaches is the handling of the nondeterminism. From the formal point of view, the non-determinism corresponds to a random equiprobable choice of an element from the set of all applicable multisets of rules (𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿)). In the case of indirect approaches, due to the iterative nature of the algorithms, it is not easy to argue that each possibility has the same probability to occur. We would state that solutions containing a smaller number of different rules have a higher chance to be selected. In the case of DND algorithm and related variants, it looks like the obtained solution tends to be an equiprobable choice. However, the corresponding articles do not give such a proof and there are some unclear points, which do not allow us to affirm this fact. Up to now, the only algorithm that is performing a truly non-deterministic choice is the one described in [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF][START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF]. However the corresponding implementation is limited to some particular derivation modes and particular classes of P systems.

AN OVERVIEW OF EXISTING FPGA SIMULATIONS

With the advent of reconfigurable hardware which realizes the idea of modifying the hardware circuits by programming, conceiving a novel circuit simulating an innovative processing paradigm is no longer an exceedingly hard task. The first attempt to use FPGA reconfigurable hardware to simulate P systems dates back to 2003 [START_REF] Petreska | A Reconfigurable Hardware Membrane System[END_REF]. Since then two simulation approaches emerged, considering regions or rules as basic processing units.

Region-based Simulations

In the region-based simulation approach rules and objects from different membranes are physically located in different places of the circuit, while those from the same membrane are physically close and well connected. The biggest problem is to ensure the correct communication of objects between membranes as this requires a global level synchronization. As advantage, the obtained system is highly scalable and robust. Below, we give two examples of region-based simulations.

In contrast, the rule-based implementation approach discussed later explicitly represents the evolution rules as processing units and multisets of objects as register arrays, while membranes and regions are represented implicitly as logical constructions existing between those processing units and data structures.

6.1.1 Petreska and Teuscher simulation. Petreska and Teuscher designed the first FPGA circuit simulating membrane computing (more precisely transitional P systems) [START_REF] Petreska | A Reconfigurable Hardware Membrane System[END_REF]. They could realize several interesting features like communication to inner membranes, priorities between rules and proposed ideas how to simulate membrane creation and dissolving mechanism in integrated circuits. This work inspired the successors to engage in this challenging and breathtaking field to advance the development of hardware simulation of P systems. In their simulation they made two strong assumptions: the application of evolution rules in each membrane is not done in a maximally parallel but in a sequential manner (but still keeping a parallelism at the system level); the nondeterministic evolution of configuration is substituted by a deterministic transition following a predetermined order. Such computational step corresponds to an ILP with the subject function as a weighted sum of variables with predefined fixed weights.

In theory, membranes are borders without internal structures and material consistence. In this implementation, the membrane structures represent regions on the circuit containing the enclosed substances, i.e., the multisets of objects, evolution rules and children membrane architectures. The objects exchange among membranes is a kind of bi-directional traversing behavior. In case of the possible objects exchange between regions, the communications are made by data buses connecting to different parts of hardware representing inner membranes. To avoid the multiple buses used to connect the parent membrane to its plural children membranes, a single bus links all the children membranes before it connects to the parent membrane. Hence, the communication is limited to parent-child membranes and there is no object exchange among children membranes or non-immediate contained membranes.

The representation of the multisets of objects is implemented by using registers. Different registers just preserve different multiplicities of objects. A register does not store the objects but only a number indicating the multiplicity of each object. The order of these registers is in accordance with the lexicographic order of the alphabet of objects. The recognition of an object is indirectly realized by examining the position of the register storing the multiplicity of this object. An evolution rule defined here is in the form of 𝑢 → 𝑣 (𝑣 1 , 𝑖𝑛 𝑖 ) (𝑣 2 , 𝑜𝑢𝑡), where 𝑣 1 is the string to be sent into lower-immediate membrane labeled i, 𝑣 2 will be sent to upper-immediate membrane. The treatment employed to deal with the formulation of evolution rules is storing the rule's left-hand side and right-hand side into different registers separately. A particular module is designed to determine whether a rule is applicable. This module compares the left-hand side of a rule u with the multiset of objects w present in the current membrane. If and only if 𝑢 ≤ 𝑤, this rule is applicable and this module will generate a signal 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 = 1. Input all the Applicable signals to an OR gate, the result of this logical gate can used as a monitor to identify whether the evolution reaches halt configuration.

The transition of configurations of P system is realized deterministically and sequentially, which is different from the general model. The consecutive transformation of configurations is regarded as the evolution process. This evolution process is decomposed into micro-steps and macro-steps. The application of rules enclosed by membranes is performed in terms of a predefined sequential order. This deterministic execution of rules is conducted in micro-steps sequentially. If a selected rule is applicable, the left-hand side of the rule u will be removed. Then the right-hand sides v, 𝑣 1 and 𝑣 2 is stored in corresponding registers. The objects from the upper immediate membrane will be preserved in another register. Although the micro-steps are carried out deterministically, they are performed simultaneously in all membranes, until there are no applicable rules. The micro-steps terminates when there are no applicable rules, i.e. the halt condition is reached. All the registers are updated in line with associated rules in macro-steps.

This implementation considered and respected the priorities of applicable rules at the beginning of each micro-step. By labeling the applicable rules with higher priorities and storing the corresponding labels, applicable rules are executed in accordance with their respective priorities. Besides, two additional features of P system, the dissolution and creation of membranes, are simulated. When a rule with membrane dissolving function is applied, its contents are owned by its upper immediate membrane, setting the membrane Enable signal of the relevant membrane to "0". However, the connections and registers defining the dissolved membrane still exist. This scheme gives rise to a disadvantage that the hardware resources cannot be released. The creation of new membranes is executed in the initialization process of the P system since all the information about new membrane is known from the specification of the system. The created membranes are inactive until membrane creating rules invoke them.

Nguyen simulation.

In this implementation, a parallel computing platform simulating membrane computing based on FPGA named Reconfig-P is developed [START_REF] Van Nguyen | An Implementation of the Parallelism, Distribution and Nondeterminism of Membrane Computing Models on Reconfigurable Hardware[END_REF][START_REF] Van Nguyen | An extensible, maintainable and elegant approach to hardware source code generation in Reconfig-P[END_REF]. Reconfig-P is fabricated on the basis of the region-oriented idea that regions work as the computational entities communicating objects through message passing. The functionality of these regions is extended by the included set of evolution rules. P Builder, the software component of Reconfig-P, specifies the P system concerned in software, converts the specification of P system written in Java to Handel-C (a hardware description language) source code. Software simulation of the circuits to be constructed is supported by P Builder to test the functionality of circuits before mapping the code to hardware circuits.

The execution of a evolution step is divided into two phases: object assignment phase and object production phase [START_REF] Van Nguyen | An Implementation of the Parallelism, Distribution and Nondeterminism of Membrane Computing Models on Reconfigurable Hardware[END_REF]. The maximal instance of each rule in a region is determined in the object assignment phase. The update of multiplicity of objects is accomplished in the object production phase. The maximal instance of the rules with higher priorities is computed before the rules with lower priorities. Note that the consumption of objects for rules with higher priorities performed during the object assignment phase to save clock cycles. It is assumed that all rules are assigned relative priorities. The priority between rules is implemented as the temporal order which should be respected by region processing units in the assignment phase. Rules with same priority are executed concurrently. The temporal order is determined at compile-time. The rules are applied according to their priorities in rounds until no rules are applicable using the indirect iterative approach. Under this circumstance, the applicability of each rule is non-stationary because of the existence of priorities. To avoid processing inapplicable rules, the applicability status of each rule is checked at the outset of the assignment phase and immediate after a applicable is applied to consume some objects.

The objects traversing behavior is the origin of communication between regions. The update of multiplicity of objects caused by rules with and without traversing behavior is completed in the object production phase. When different region processing units update the multiplicity value of the same object at the same time, a conflict occurs. To handle this conflict, in [START_REF] Van Nguyen | An Implementation of the Parallelism, Distribution and Nondeterminism of Membrane Computing Models on Reconfigurable Hardware[END_REF] two solution strategies, the space-oriented strategy and the time-oriented strategy are proposed. Table 1 and2 summarize the different strategies of two resource conflict resolutions and their modifications in the rule-based and region-based design [START_REF] Van Nguyen | Balancing Performance, Flexibility, and Scalability in a Parallel Computing Platform for Membrane Computing Applications[END_REF][START_REF] Van Nguyen | An Implementation of Membrane Computing Using Reconfigurable Hardware[END_REF]. In order to simplify the exposition of processes of rule-based and region-based implementations of P systems, tables are designed to delineate the relevant details, which will be given below. The number of delays is equal to 𝑠. Space-oriented strategy (a) construct the conflict matrix as in the Time-oriented strategy. (b) the register storing the object accessed by multiple rules concurrently is replicated 𝑠 times. These copy registers are assigned to each conflicting rules to write. After the updating process for all the copy registers, the corresponding values are joined to the original register.

Table 2. The differences of the conflict resolutions adopted in two design modes Item Time-oriented strategy Space-oriented strategy

Rule-oriented design

The interleaving operation can be determined at compile-time and it can be hard-coded into the HDL source.

Need a multiset replication coordinator to coordinate the multiplicities stored in the copy registers.

Region-oriented design The objects received from other regions are regarded as external objects, otherwise the objects are internal objects. The interleaving merely caused by the production of internal objects can be identified at compile-time. To retain the independence of region processing units, the interleaving induced completely or partially by the receipt of external objects can only be calculated at runtime.

The role of the multiset replication coordinator in ruleoriented design is played by the region processing units. The existing register storing the multiplicity received from the associated communication channels in the considered region can be assigned to those processing units which sending objects to the considered region to write the new values of of the objects competed by multiple rules.

The extensibility of the region-based design is the consequence of the representation of membranes as processing units interacting with two region processing units corresponding to inner and outer regions. This allows to achieve a strong separation of the processing logic inside different membranes and the independence of the communication. Thus, adding additional elements to the system does not lead to the redesign of the remaining part of the system.

Rule-based simulations

Rule-based approaches consider evolution rules as processing units performing the update of multiplicities and of the membrane structure.

Nguyen simulation .

Every rule in all regions of the P system is represented as a processing unit synchronized by a global clock that implements the parallel processing. However, a processing unit does not correspond a concrete hardware component but rather to a potential infinite while loop which contains the code related to the rule application using Handel-C. The information associated to execution and synchronization is contained in processing units as well. Each rule processing unit in a region is linked to the array of registers containing multisets of objects. The membrane inclusion relationships can be described with the connections between processing units and arrays. Generally speaking, a rule processing unit in a region is linked to the objects array located in the same region with the rule processing unit. If there are objects traversing regions which imply the containment, connecting the rule processing unit to the object array to which the rule will send objects contained in different region. This procedure permits to represent the membrane inclusion. The rule execution is split into preparation phase and updating phase. We compare the two designs (the rule-based and the region-based) of Nguyen's implementation in Table 3. Regions are realized in hardware implicitly by the content they included. For the containment relationships including the traversing of objects between regions, they are implemented by imparting the corresponding rules with 'in' or 'out' target directives the abilities to access the multiset of objects in the destination region of the objects traversed.

Regions are represented as parallel processing units. The traversing of objects among regions is realized as message passing through channels connecting different region processing units.

Multiset of objects An array of registers contains each type of object in the alphabet for a region

The same strategy adopted as the rule-oriented design.

Evolution rule Potentially infinite while loop in which contain procedures representing the operations of the relevant evolution rules in HDL language.

Implicitly expounded through integrating them into the region processing units.

Operation process Preparation phase and updating phase. Object assignment phase and production phase.

Synchronization

An array of registers composed of three 1-bit registers is associated to each rule processing unit. The values in the array indicate whether the rule processing unit has complete the preparation and updating phase, or is applicable. Compute the logic AND or OR values of all the values stored in the 1-bit register that indicates the same status of the processing unit and store the three results into three sentinel registers. The coordinating processing unit read the sentinel values to synchronize the whole procedures.

For the synchronization of object assignment phase, it is realized when all the region processing units communicate with each other on channels at the beginning of object production phase. For the object production phase, a region execution coordinator connecting to each region processing unit via dedicated channels is designed to perform the synchronization of operations.

After the region execution coordinator received the signals denoting the completion of all the operations of every region processing unit, the current transition is done and the next one is carried out.

Verlan and Quiros simulation .

The target model is a static P system. The system is considered flattened, so only one skin membrane is present. A special strategy was elaborated in order to not compute the complete solution, i.e. 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿) but the cardinality of its elements. Then a random value between 1 and this cardinality is taken. Finally, this number is decoded to the corresponding solution [START_REF] Quiros | Fast Hardware Implementations of Static P Systems[END_REF][START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF] Devising an algorithm which carries out the computation of the cardinality and of all elements of the solution set in constant time on FPGA is the key issue of the approach. A remarkable characteristic of FPGA is that the time consumed for executions of functions that do not exceed the cycle of the global clock is done in one cycle of FPGA, hence in constant time. The computation of the cardinality and the decoding of solutions are accomplished by two functions hardwired into the circuit: 𝑁 𝐵𝑉 𝑎𝑟𝑖𝑎𝑛𝑡𝑠 (Π, 𝐶, 𝛿), which gives the cardinality of 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿), and 𝑉 𝑎𝑟𝑖𝑎𝑛𝑡 (𝑛, Π, 𝐶, 𝛿), which returns the 𝑛-th element of 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝛿).

A concept named rules' dependency graph is introduced to compute the two functions above. It is a bipartite graph that contains as nodes rules and objects and there is an edge between two nodes if a rule contains the corresponding objects in its left-hand-side. The picture below depicts the rules' dependency graph for rules 𝑟 1 : 𝑎𝑏 → 𝑢 and 𝑟 2 :

𝑏𝑐 → 𝑣. 𝑟 1 𝑟 2 𝑎 𝑏 𝑐 d d d d
Assume that the derivation mode is maximal parallelism (max). Suppose that 𝑁 𝑎 , 𝑁 𝑏 and 𝑁 𝑐 represent the number of objects 𝑎, 𝑏 and 𝑐 in 𝐶. Let 𝑁 1 = 𝑚𝑖𝑛(𝑁 𝑎 , 𝑁 𝑏 ), 𝑁 2 = 𝑚𝑖𝑛(𝑁 𝑏 , 𝑁 𝑐 ), 𝑁 = 𝑚𝑖𝑛(𝑁 1 , 𝑁 2 ), 𝑘 𝑖 = 𝑁 𝑖 ⊖ 𝑁 ,1 ≤ 𝑖 ≤ 2, where ⊖ denotes the positive subtraction. Let also 𝑝, 𝑞 = 0, 1, 2, . . . , 𝑁 . From the dependency graph we can deduce the following:

𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) = 𝑝+𝑞=𝑁 𝑟 𝑝+𝑘 1 1 𝑟 𝑞+𝑘 2 2 𝑁 𝐵𝑉 𝑎𝑟𝑖𝑎𝑛𝑡𝑠 (Π, 𝐶, 𝑚𝑎𝑥) = 𝑁 + 1 𝑉 𝑎𝑟𝑖𝑎𝑛𝑡 (𝑛, Π, 𝐶, 𝑚𝑎𝑥) = 𝑟 𝑁 -𝑛+1+𝑘 1 1 𝑟 𝑛-1+𝑘 2 2 .
Example 6.1. Consider a configuration where 𝑁 𝑎 = 5, 𝑁 𝑏 = 5 and 𝑁 𝑐 = 3. It can be easily verified that

𝑁 1 = 𝑚𝑖𝑛(5, 5) = 5, 𝑁 2 = 𝑚𝑖𝑛(5, 3) = 3, 𝑁 = 𝑚𝑖𝑛(5, 3) = 3, 𝑘 1 = 𝑁 1 ⊖ 𝑁 = 5 -3 = 2, 𝑘 2 = 𝑁 2 ⊖ 𝑁 = 3 -3 = 0.
Hence, we can enumerate the elements of 𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) as bellow:

𝐴𝑝𝑝𝑙𝑖𝑐𝑎𝑏𝑙𝑒 (Π, 𝐶, 𝑚𝑎𝑥) 1 = {𝑟 3+2 1 𝑟 0+0 2 , 𝑟 2+2 1 𝑟 1+0 2 , 𝑟 1+2 1 𝑟 2+0 2 , 𝑟 0+2 1 𝑟 3+0 2 } = {𝑟 5 1 , 𝑟 4 1 𝑟 2 , 𝑟 3 1 𝑟 2 2 , 𝑟 2 1 𝑟 3 2 }.
The same result can be easily obtained by using formal power series associated to context-free languages. In this case, any maximal rule combination is a part of the language

𝐿 𝑁 = {𝑟 𝑝 1 𝑟 𝑞 2 | 𝑝 +𝑞 = 𝑁 }.
It is quite easy to observe that the number of words of length 𝑁 in 𝐿 𝑁 is exactly the same as the number of words of same length in the language 𝐿 = {𝑟 * 1 𝑟 * 2 }. This last language is regular and its generating function is 𝑞 0 (𝑥) = 1/(1 -𝑥) 2 . The 𝑛 -𝑡ℎ coefficient of the expansion of 𝑞 0 (𝑥) is equal to 𝑛 + 1 ([𝑥 𝑛 ]𝑞 0 = 𝑛 + 1), which immediately gives 𝑁 𝐵𝑉 𝑎𝑟𝑖𝑎𝑛𝑡𝑠 (Π, 𝐶, 𝑚𝑎𝑥) = 𝑁 + 1. The 𝑉 𝑎𝑟𝑖𝑎𝑛𝑡 (𝑛, Π, 𝑐, 𝑚𝑎𝑥) is computed using an algorithm that performs a weighted breadth-first search of the decomposition of 𝑛 with respect to the number of variants found on each branch of the execution of automaton for 𝐿.

Such a process can be easily repeated for any regular language, yielding a constant time simulation of a computational step. The reason for such performance is that any generating function is equivalent to a recurrence relation and such relations can be computed in one synchronous time unit using asynchronous operations. The described algorithm functions for any P system where the rule choice can be expressed as words of certain length in a regular language. The corresponding class is quite large (containing even computationally complete models), thus allowing an extremely fast execution. Examples from [START_REF] Verlan | Fast Hardware Implementations of P Systems[END_REF] exhibited a speed-up of order 10 5 .

Another important point is that this approach allows to handle the non-determinism in a natural way, by performing a uniform random choice between all possible rule applications at each step.

Technically, the implementation represents only objects by registers and rules by layered logic. Each rule implementation is modularized and contains an own copy of processing instructions needed to compute the two above functions, based on asynchronous operations. Consequently, five clock cycles are required to compute the 𝑁 𝐵𝑉 𝑎𝑟𝑖𝑎𝑛𝑡𝑠 (Π, 𝐶, 𝑚𝑎𝑥), 𝑉 𝑎𝑟𝑖𝑎𝑛𝑡 (𝑛, Π, 𝐶, 𝑚𝑎𝑥) and to apply the corresponding rules. The entire process of the implementation is split into several consecutive stages, which take charge of different operations associated to phases of evolutions of configurations.

Persistence stage stores the states that the hardware system goes through. A independent stage computes the maximal instance of each rule by means of the dividing operation and MIN logic operation. Assignment stage in charge of selecting a rule to be applied non-deterministically, and determines its instances. Updating stage is responsible for updating the current configuration with the values from the previous stage. During Halting stage, the system inspects whether the halting condition is reached, and once reached, stops the system.

According to the elaborations stated above, the comparisons of the FPGA implementations of P systems is listed in table 4. The details of each design are summarized. We can check the characteristics of the three main prototypes. Subjects Cell-like P systems with following characteristics: membrane dissolution and creation, objects exchange between upper-and lowerimmediate membranes, cooperative P systems with priorities. The range (non-continuous) for object number is 6 to 12 and for membrane is 10 to 20.

For rule-oriented design, the subjects are cell-like P systems cascaded in vertical, horizontal, vertical and horizontal structures. For regionoriented design, the objects are celllike P systems containing hierarchical regions and tissue-like connected regions. The rule range is [START_REF] Jambhlekar | Parallel implementation of MOPSO on GPU using OpenCL and CUDA[END_REF][START_REF] Kaliszewski | Multiple Criteria Decision Making by Multiobjective Optimization-A Toolbox[END_REF], [START_REF]IRDS points to chip stacks and new architectures[END_REF][START_REF] Díaz-Pernil | Membrane computing and image processing: a short survey[END_REF] for regions, [START_REF] Agrigoroaiei | Evolving by Maximizing the Number of Rules: Complexity Study[END_REF]200] for objects. The extent of the inter-region communication is [START_REF] Pérez-Hurtado | P-Lingua in two steps: flexibility and efficiency[END_REF]319]. The subject P system is simplified according to the multiset rewriting point of view, which it has a skin membrane, no inner regions. The 4 subjects differ in the rule dependencies which form chains: circular, 2circular, linear, opposite. The object number range is [START_REF] Jambhlekar | Parallel implementation of MOPSO on GPU using OpenCL and CUDA[END_REF]200].

Experiment results

The hardware consumption ranges from 4.2% to 33% (CLB). The extent of clock rate is 27 to 198 MHz.

For rule-oriented design, the number of rules applied per second ranges from 2.7×10 5 to 1000×10 5 . The hardware consumption extent is 1.55% to 21.43% (LUTs). For region-oriented design, the hardware usage ranges 1.82% from 16.79%. The clock rate fluctuates from 52.63 MHz to 81.77 MHz. The biggest size which can be executed is a P system with 550 rules, 1280 communication channels, 1100 objects conflicts.

The hardware consumption ranges from nearly 1% to 48% (LUTs) or nearly 1.1% to 11% (slice). The period needed to perform a computation step fluctuates from 5.46 ns to 9.14 ns. The highest frequency exceeds 100 MHz, permitting 2 × 10 7 computational steps per second. The run-time for each experiment subject ranges from 3.017 × 10 -5 𝑠 to 4.174 × 10 -5 𝑠.

Parallelism System-level parallelism Region-level and system-level parallelism System-level (there is only a skin membrane, so it is also region-level parallelism).

Non-determinism

No non-determinism DND algorithm True non-determinism.

AN OVERVIEW OF EXISTING CUDA SIMULATIONS

In [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF], it is concluded that the GPU is a suitable platform to accelerate the simulation of P systems because of the following features:

• Good performance: for example, the NVIDIA Tesla K40 delivers 1.43 TeraFLOPS doubleprecision peak floating point performance, 4.29 TeraFLOPS of single-precision, and 288 GBytes/s of global memory bandwidth; • An efficiently synchronized platform: GPUs implement a shared memory system, avoiding communication overload; • A medium scalability degree: the amount of resources depend on the GPU model, e.g. a K40

includes 2880 cores and 12 GBytes of memory. If the resources of a GPU is not enough, there are more scalable solutions such as multi-GPU systems, but they then require communication among nodes; • Low-medium flexibility: although CUDA programming is based on C++, and hence programmers are free to use the same data structures than in CPU, both the algorithm and the data structure have to be adapted for best performance on GPUs.

Simulation algorithms implemented on CUDA have a common structure, in which for each simulated computation step, first the rules are selected (obtaining a multiset of rules, while consuming the left-hand sides), and secondly the rules are executed (based on the obtained multiset of rules, Membranes (regions) and their containment Implicitly represented by the contents enclosed by the membranes. The objects exchange is interpreted as transferring objects with communication buses by which connect the origin region to the destination region.

For the rule-oriented design, the treatment is similar with Petreska's. For the region-oriented design, regions are represented as parallel processing units. The objects exchange among regions is realized as message passing through channels connecting different region processing units.

According to the multiset rewriting system framework, the topology structures of membranes are not important. What concerned is the rule dependency graph.

Multiset of objects Store the multiplicity value of each type of object in different registers whose positions indicate the type.

The same method. The same method.

Evolution rule Store the left-hand-side and righthand-side of a rule in different registers.

In rule-oriented design, they are characterized as potentially infinite while loop in which contain procedures representing the operations of the relevant evolution rules in HDL. In region-oriented design, they are implicitly expounded through integrating them into the region processing units.

The logic of rules is distributed along the the hardware components. There is no explicitly correspondence between rules and hardware components.

Operation Process Micro-step: only one applicable rule is applied with respect to the instance number in a region. Macro-step: execute a micro-step concurrently in every region.

In rule-oriented design, perform preparation phase and updating phase. In region-oriented design,perform object assignment phase and production phase.

1. Persistence stage 2. Independent stage 3. Assignment stage 4. Application stage 5. Updating stage 6. Halting stage.

Extensibility

Membrane mediated features cannot be added in since membranes are implicitly represented by their contents.

Because of the region-oriented design, Membrane mediated features such as symport and antiport functions can be extended.

This implementation is designed only for P systems whose applicable multisets of rules can be represented as regular language. So its extensibility is limited.

Scalability

With the increase of the rules, the hardware usage rise approximately proportional. The clock rates decline a small amount. The membrane creation ability cripples clock rate significantly.

The hardware consumption scales linearly with respect to the size of the P system executed in both rule-oriented and region-oriented design. The performances grow linearly when the number of rules increase.

The hardware usage is scalable, but the rate of increase in the performance is not always linear for different systems with distinct structures.

Contributions Membrane creation and dissolution functionality

Two kinds of methodology for P system characterization and resource conflict resolution, DND algorithm Put up with a new methodology with absolute equiprobability to implement non-determinism.

Drawback

Partially parallelism, no nondeterminism

The equiprobability of DND algorithm has not been proven theoretically.

The types of P systems which can be implemented are confined to those whose rules can be represented as regular language.

and generating the right-hand side. This strategy is necessary in order to synchronize which rules get executed in the corresponding computation step. For some models, for example for those ad-hoc simulators, the selection step is done in micro-stages.

In the following subsections, the existing CUDA simulations are summarized, by organizing into P system models.

Cell-like P systems

The first test of concept for simulating P systems on GPUs was applied to P systems with active membranes using CUDA [START_REF] José | Simulation of P systems with active membranes on CUDA[END_REF]. This simulator performs only one computation out of the whole tree, in order to avoid non-determinism, by requiring the confluence property to the simulated P systems. Bearing this in mind, the "lowest-cost computation path" is selected: the one in which least membranes and communication are required. This is achieved by giving preferences to rules that lead to least membranes (e.g. dissolution over division rules). The simulation algorithm is composed of two main stages: selection and execution of rules. Selection is where the semantics of the model is actually simulated. Rules are chosen by following the defined constrains, altogether with a number of applications. The result of this stage is used for the next one, which is the execution of the rules; that is, updating the P system configuration. This two-staged strategy allows to synchronize the application of rules within and among membranes.

Both P systems and GPUs have a double-parallel nature [START_REF] José | Simulation of P systems with active membranes on CUDA[END_REF], and this is harnessed for implementing a mapping: (elementary) membranes are assigned to thread blocks, and a subset of rules to threads. Each thread is in charge of selecting rules for a portion of the defined objects in the alphabet. Note that this is enough given that in P systems with active membranes, rules have no cooperation. However, this mapping of parallelism is naive, since it assumes that all the objects in the alphabet can be present within each membrane. This require allocating memory space and assigning resources (threads) to all of them. This, in fact, does not take place in the majority of P systems to be simulated, but turns out to be the smallest worst case to handle. Thus, the performance of the simulator completely depends on the P system being simulated, and drops as long as the variety of different objects appearing in membranes decreases.

Non-determinism is handled by imposing the simulated P systems to be confluent; that is, all computations halt and they generate the same result [START_REF] José | Simulation of P systems with active membranes on CUDA[END_REF]. This way, the simulator can choose any path in the computation tree, since the aim is to find a halting configuration. The GPU simulator take advantage of this by preferring to select evolution rules rather than division or dissolution.

The performance of the simulator was analyzed on a GPU Tesla C1060 (240 cores, 4GB memory) by using two benchmarks [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF]: a simple test P system designed in a convenient manner (up to 7x of speedup), and a family of P systems designed to solve different instances of the SAT problem (1.67x of speedup).

Improvements to this design have followed [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF], reporting up to 38x of acceleration when taking advantage of shared memory and data transfer minimization. By constructing a dependency graph, the set of rules of the input model are arranged so that those having common objects in the lefthand side and in the right-hand side (respectively) are more likely to be in a node. This reduces communication given that thread blocks are assigned to nodes.

Another approach was to implement ad-hoc simulators for a specific family of P system with active membranes allowing to solve SAT in linear time [START_REF] José | The GPU on the simulation of cellular computing models[END_REF]. In this way, the worst case assumed before (all objects defined can appear in every membrane) is further reduced by analyzing the upper bound to the number of existing objects in membranes. In this way, the work done by threads is maximized, and the design remains similar: each elementary membrane is implemented by a thread block, and each object of the input multiset is implemented by a thread. The experiments carried out on a NVIDIA Tesla C1060 GPU reported up to 63x of speedup. Further developments took place focusing on this family of simulators, with the aim at being better tailored to newer GPU architectures, and also to enable multi-GPU systems and supercomputers [START_REF] José | The GPU on the simulation of cellular computing models[END_REF].

A related work was to explore which P system ingredients are better suited to be handled by GPUs. To this aim, another solution to SAT based on a family of tissue P systems having the operation of cell division was simulated [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF]. The design is similar to the previous case: each cell is simulated by a thread block, however the constructed solution required a higher number of objects to be placed inside each cell. At the same time, this simulator does not need to store nor handle charges associated to membranes. Experiments on a NVIDIA Tesla C1060 GPU leaded to speedup by 10x. This showed that using charges associated to membranes helped to save instantiation of objects, and so, they entail a lightweight ingredient to be processed by threads.

Population Dynamics P systems

Population Dynamics P (PDP) systems is a multi-environmental model successfully used for modeling real ecosystems. Thus, their efficient simulation was critical for experimental validation.

Simulators for PDP systems typically run several simulations in order to extract statistical information from the models.

A CUDA simulator for PDP systems was presented in [START_REF] Miguel | Parallel simulation of Population Dynamics P systems: updates and roadmap[END_REF]. The selected simulation algorithm was the DCBA [START_REF] Miguel | DCBA: Simulating Population Dynamics P Systems with Proportional Object Distribution[END_REF], since it provides better accuracy in the simulation results. The selection of rules in DCBA consists of three phases: phase 1 (distribution of objects), phase 2 (maximality) and phase 3 (probability). This algorithm uses a distribution table in order to distribute the objects in a proportional way between competing rules, i.e. with overlapping left-hand sides. Moreover, rules are grouped into rule blocks, when having the same left-hand side.

It can be noticed that the approach of DCBA is to handle non-determinism by first using an uniform distribution of objects to competing rules. This requires an extra phase (2) to avoid rounding errors, in this way a random order is employed to control the remaining rules and which one takes all objects. This phase is based on the procedure of DND.

The CUDA simulator for PDP systems [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF][START_REF] Miguel | Parallel simulation of Population Dynamics P systems: updates and roadmap[END_REF][START_REF] Miguel | Adaptative parallel simulators for bioinspired computing models[END_REF] uses the following design: the contents of environments and rule simulations are fractioned and distributed all-over thread blocks, while individual rules are fractioned over threads. Phases 1, 3 and 4 are efficiently managed by the GPU, but Phase 2 can become a bottleneck. For phase 3, a random binomial variate generation library was developed, so that binomial and multinominal distributions were supported for random number generation. In the benchmark using a set of randomly generated PDP systems (having no biological meaning), speedups of up to 7x were achieved on a Tesla C1060 with respect to a multi-core version. The simulator was validated and tested by using a known ecosystem model of the Bearded Vulture in the Catalan Pyrenees, leading to speedups of up to 4.9x with a C1060, and 18.1x using a Tesla K40 GPU (2880 cores). Finally, the simulator was extended in [START_REF] Miguel | Adaptative parallel simulators for bioinspired computing models[END_REF] in such a way that it introduces high-level information provided by the model designer into the code. This is accomplished through a new syntactical ingredient in P-Lingua 5 [START_REF] Pérez-Hurtado | P-Lingua in two steps: flexibility and efficiency[END_REF], called feature. This ingredient was employed to define the usual algorithmic scheme based on rule modules when designing biological models. This new version is called adaptative simulator, and helped to obtain an extra speedup of 2.5x on a Tesla K40, and 1.8x on a Tesla P100 GPU (3584 cores).

Spiking Neural P systems

Parallel simulation of Spiking Neural P (SNP) systems has been based on a matrix representation so far [START_REF] Zeng | Matrix Representation of Spiking Neural P Systems[END_REF]. The simulation algorithm uses the following vectors and matrices:

• Spiking transition matrix: stores information about rules, and is employed for computing transitions. It assigns a row per rule and a column per neuron. The values coded correspond to the left-hand side of rules (as negative numbers) and to the right-hand sides (as positive numbers). • Spiking vector: defines a selection of rules to be fired in a transition step, using a position per rule. Given the non-deterministic nature of SNP systems, there are more than one valid spiking vector for a given configuration. • Configuration vector: defines the number of spikes per neuron, that is, the configuration in a given time.

The algorithm is then used to compute the next configuration from a given one, by performing only vector-matrix operations. These operations are offloaded to the GPU, which is optimized to handle matrices [START_REF] Francis | Improving GPU simulations of Spiking Neural P systems[END_REF]. CuSNP is a simulator for SNP systems which is written in Python and the CUDA kernels were launched by using the binding library PyCUDA. For the first approach, SNP systems without delays were simulated by covering each computation path sequentially, leading to speedups of up to 2.31x [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF] [START_REF] Francis | Improving GPU simulations of Spiking Neural P systems[END_REF].

Further extensions have followed, enabling delays, support of more types of regular expressions and input of P-Lingua files [START_REF] Jym | CuSNP: Spiking Neural P Systems Simulators in CUDA[END_REF]. The support of delays is done by an extension of the matrix representation, introducing vectors to control when rules fire after the delay is met, and if neurons are open. These leaded to up 50x of acceleration on a GTX750 GPU [START_REF] Jym | CuSNP: Spiking Neural P Systems Simulators in CUDA[END_REF]. Moreover, non-deterministic SNP systems with delays are supported by generating all computation paths sequentially and using GPU to speedup the simulation [START_REF] Jym | Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations[END_REF] (with up to 2x on a GTX1070 for a non-uniform solution to Subset Sum).

Furthermore, the simulation of Fuzzy Reasoning Spiking Neural (FRSN) P systems on the GPU was explored [START_REF] Luis | Simulating FRSN P Systems with Real Numbers in P-Lingua on sequential and CUDA platforms[END_REF]. FRSNP systems allows modeling fuzzy diagnosis knowledge and reasoning for fault diagnosis applications. The simulation algorithm is also based on a matrix representation and vector-matrix operations. The employed simulation framework was pLinguaCore, so the CUDA kernels were launched by using the binding library JCUDA.

Other models

Enzymatic Numerical P systems has been employed for modeling robot controllers and path planning, being significant for the Artificial Intelligence. A CUDA simulator was developed [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF], where the selection of applicable programs was first applied, second, the calculation of production functions, and third, distribution of production function results according to repartition protocols. Production functions are computed using a recursive solution. Simulators were implemented in Java (inside pLinguaCore) and C programming languages as standalone tools. On a GeForce GTX 460M, the achieved speedup was of up to 11x. Moreover, a model for RRT and RRT* algorithms for robotic motion planning using ENPS was also simulated on CUDA [START_REF] Pérez-Hurtado | A membrane parallel rapidly-exploring random tree algorithm for robotic motion planning[END_REF] 24x faster than in multicore CPU using an RTX2080.

Evolution-Communication P systems with Energy (ECPE) [START_REF] Miguel | Simulating P Systems on GPU Devices: A Survey[END_REF] were also target for a CUDA simulation. The employed simulation algorithm used a matrix representation and linear-algebra based algorithm, similarly as for the spiking neural P systems simulator: a configuration vector, a trigger matrix, an application vector and a transition vector.

The summary of CUDA-GPU implementations is concluded in Table 6 and7, providing an easy checking for the CUDA-based simulations for academic communities. For the sake of presenting a overview contrast of FPGA and CUDA-GPU implementation, we summarize the quantitative and qualitative attributes and conclude their methodologies in the Table 8, 9 and 10. This overall contrast present a recapitalized conclusion about hardware implementations of P system. Programming with HDL is difficult contrasting with its software counterpart, and the synchronization of intricate circuits. It usually cost 8-10 months to be proficient with HDL programming and to comprehend the parallel architecture.

Work assignation for a parallel design can be the most difficult task when developing GPU simulators, since the restrictions of GPU computing meets with P system semantics. Need to think on GPU architecture for high performance: memory accesses, memory transfers, thread synchronization, etc. Application

The prototypes are developed but no applications for concrete problems to date.

several application for biological simulations have to executed. Solutions to computationally-hard problems.

Table 10. The conclusion of the implementations for FPGA and CUDA-GPU

Item FPGA CUDA-GPU

Advantages

The computation speed is fast and the hardware framework is adjustable. The compromises between implementing the P system model and constructing the most proximate hardware are relatively ideal.

Relatively convenient way to implement P systems for the established parallel framework in the GPU. The programming is similar with the classic software developing process.

Disadvantages

Building the parallel architecture from scratch is a laborious and challenging adventure Slower than FPGA. Sometimes big concessions are indispensable as the architectures are unchangeable.

Contributions

Introducing the reconfigurable hardware to develop the real parallel architectures which can exploit the maximal parallelism of the P systems substantially.

The tremendous potential of the emergent universal computing GPU is concentrated, which the readymade parallel framework is off-the-shelf, to establish a P system computing platform conveniently.

Conclusions

We can develop sophisticated P systems hardware circuits which implement the target models as approximate as possible on FPGA devices. The expense we should pay is also considerable high taking the painstaking effort into account.

CUDA-GPU provides a relatively comfortable and alternative choice to implement P systems, which the simulation results are acceptable, nevertheless at a price of the modification of the source models. GPUs worth when simulating large P system models.

OTHER APPROACHES

Besides FPGA-based and CUDA-enabled GPU-based hardware platforms developed for the implementation of P systems, there are several attempts to implement P systems on micro-controllers. Serial algorithms and their hardware circuits designs in terms of the exhaustive investigation line focusing on implementing the transition of configurations of P systems were developed. The corresponding research does not precise concrete hardware devices, it is just designing the circuits aiming at simulating certain operations of particular P systems with registers, logical gates, magnitude comparators, and data buses.

In [START_REF] Fernández | A Hardware Circuit for Selecting Active Rules in Transition P Systems[END_REF], a digital circuit is presented to select active rules in the current configuration. Each evolution rule is represented by 2 hardware registers. The first register characterizes the left-hand side (antecedent) of a rule, and the other specifies the right-hand side of the rule, which also determines whether the rewritten objects go out from the current membrane or stay where they are, or go into inner membranes. By comparing the left-hand side of each rule with the multiset of objects in a region, the applicability of every rule can be determined. Next, an algorithm computes the number of applications of active rules given the multiset of objects and evolution rules [START_REF] Martínez | Hardware Implementation of a Bounded Algorithm for Application of Rules in a Transition P-system[END_REF]. The corresponding circuit is composed of logical gates, registers, multiplexers and sequential elements. The computation process is bounded. In [START_REF] Gutiérrez | Design of a Hardware Architecture Based on Microcontrollers for the Implementation of Membrane Systems[END_REF], a P system circuit is constructed by means of PIC16F88 micro-processor enriched with the storage component 24LC1025, connected by an I2C bus. Thus, the shortcoming of insufficient storage capacity of the micro-processor is solved by the introduced external memory. The flexibility of the circuit is acceptable as the modification of the structure is not necessary.

Continuing the research from [START_REF] Martínez | Hardware Implementation of a Bounded Algorithm for Application of Rules in a Transition P-system[END_REF], the paper [START_REF] Martínez | HW Implementation of a Optimized Algorithm for the Application of Active Rules in a Transition P-system[END_REF] presents an improved algorithm and the corresponding circuit calculating the application times of active rules. The computing process can be completed in minor steps and the theoretical performance is optimized. In [START_REF] Alonso | A Circuit Implementing Massive Parallelism in Transition P Systems[END_REF][START_REF] Alonso | Main modules design for a HW implementation of massive parallelism in transition P-systems[END_REF], a draft of a circuit implementing the inherent parallelism of P system is presented. Towards the parallelism, the rule treatment is similar with the region-based solution defined in Section 6, namely, what applied is a multiset of rules instead of a single rule. An operating environment is elaborated in [START_REF] Sandra | Hardware Implementation of P Systems Using Microcontrollers. An Operating Environment for Implementing a Partially Parallel Distributed Architecture[END_REF] which performs the automatic transformation of tasks involved in the hardware simulation of P systems, including loading, execution and interpretation, into a distributed framework constructed on micro-controllers. The execution results of the circuit, which are in the form of binary data, can be interpreted in a transparent manner. What should be emphasized here is that all the hardware circuits introduced are just on the blueprints which were never implemented in practice. They remain mostly theoretical and the actual functionality and performances are unknown, unlike for FPGA-based and CUDA-based hardware implementations/simulations which are carried out practically. In [START_REF] Gutiérrez | Suitability of using microcontrollers in implementing new P-system communications architectures[END_REF], micro-controllers are also chosen as target hardware to implement communication architectures of P systems. A digital circuit carrying out massive parallelism in transition P Systems is established in [START_REF] Alonso | A Circuit Implementing Massive Parallelism in Transition P Systems[END_REF].

As a new attempt for implementing neural P systems on different hardware, DRAM-based CMOS circuits are adopted to construct elementary spiking neural P systems. We do not carry out an in-depth discussion about this topic and refer to [START_REF] Xu | The stochastic loss of spikes in spiking neural P systems: Design and implementation of reliable arithmetic circuits[END_REF] for more details.

CONCLUSIONS

In this article we gave an overview of different existing hardware implementations of P systems. As expected the FPGA-based approach is very promising, yielding truly parallel implementations achieving a speed-up of order 10 5 . However, it requires a considerable amount of work as the corresponding hardware architecture should be created from scratch. Moreover, the hardware description languages like VHDL are very low-level and programming in such languages tends to be extremely tedious. Also, as it was shown from the existing implementations, the most promising results are obtained using rule-based architectures, because of the smaller communication cost.

By contrast, CUDA-based approach features a powerful unique hardware platform that can be programmed using standard C/C++ language. While it gives less freedom and accuracy, this approach gives a good trade-off between complexity, scalability and maintenance. The achieved speed-ups are relatively small (between 4 and 50 times), however, due to the programming simplicity, more types of P systems were implemented.

The central problem in both approaches -the object distribution problem -was tackled from different points of view, the most fruitful attempts being the variations of the DND algorithm and direct approaches using mathematical properties of the dependency relation between rules. Future research can be done in this direction by providing new classes of P systems suitable for the precomputation of possible rule applications. Another research direction related to software implementation of P systems is given by the construction from Section 4, which features a novel reduction of the object distribution problem to ILP. This reduction handles well the maximal parallelism and allows to define criteria for the choice of the solution. Then, it would be possible to design a software framework that would use existing solvers and algorithms to quickly obtain the desired solution, including those optimized for parallel hardware [START_REF] Jambhlekar | Parallel implementation of MOPSO on GPU using OpenCL and CUDA[END_REF][START_REF] Daniele | A Fine-Grained CUDA Implementation of the Multiobjective Evolutionary Approach NSGA-II: Potential Impact for Computational and Systems Biology Applications[END_REF][START_REF] Talbi | Parallel Approaches for Multiobjective Optimization[END_REF][START_REF] Talbi | A unified view of parallel multi-objective evolutionary algorithms[END_REF][START_REF] Tsotskas | The Design and Implementation of a GPUenabled Multi-objective Tabu-search Intended for Real World and High-dimensional Applications[END_REF][START_REF] Zhu | DEMCMC-GPU: An Efficient Multi-Objective Optimization Method with GPU Acceleration on the Fermi Architecture[END_REF].

Finally, as a future research interest we mention a robust implementation of different features of P systems like membrane creation/dissolution, as well as non-classical variants of P systems like numerical P systems. The motivation for such research is that the corresponding models feature numerous applications, so their scalable implementation would immediately allow to test the acceleration of corresponding algorithms and their practical usage.

Fig. 1 .

 1 Fig. 1. Feasible and Pareto-optimal solutions from Example 4.2.

Example 4 . 3 .

 43 Consider the system Π = (𝑂, 𝑤 1 , 𝑅), with 𝑂 = {𝑎, 𝑏, 𝑐} and 𝑅 = {𝑟 1 : (1, 𝑎𝑏) → (1, 𝑎𝑏𝑐); 𝑟 2 : (1, 𝑏𝑏𝑐) → (1, 𝑎𝑏𝑏)}. Consider the configuration 𝐶 = (1, 𝑎 5 𝑏 6 𝑐 3 ). Then, the constructed IMCLOP corresponds to the one given in Example 4.2

Example 4 . 4 .

 44 Consider the system Π = (𝑂, 𝑤 1 , 𝑅), with 𝑂 = {𝑎, 𝑏, 𝑐} and 𝑅 = {𝑟 1 : (1, 𝑎𝑏𝑐) → (1, 𝑎𝑏); 𝑟 2 : (1, 𝑎) → (1, 𝑏𝑏); 𝑟 3 : (1, 𝑏) → (1, 𝑐𝑏)}. Consider the configuration 𝐶 = (1, 𝑎 2 𝑏 3 𝑐 2

Table 1 .

 1 The comparison of the time-oriented and space-oriented conflict resolutionStrategyResource conflict resolution Time-oriented strategy (a) construct a conflict matrix in which each row is a quadruple (𝑝, 𝑞, 𝑟, 𝑠). 𝑝 is the object competed by multiple rules. 𝑞 is the region where 𝑝 is produced or consumed. 𝑟 is the set of the conflicting rules, 𝑠 is the size of set 𝑟 . (b) insert delay statement among conflict rules such that the updating operations of multiplicity of 𝑝 can be executed in distinct clock cycles.

Table 3 .

 3 The comparison of the rule-based and region-based design of Nguyen's implementation

	Object	Rule-based Design	Region-based Design
	Region and their contain-
	ment relationships

Table 4 .

 4 The quantitative attributes of FPGA implementation

	Item		Biljana Petreska	Van Nguyen	Juan Quiros
	Period		2003	2007-2010	2012-2015
	Institute		Swiss Federal Institute of Technology	University of South Australia	University of Seville
			Lausanne (EPFL)
	Target FPGA	Xilinx Virtex-II Pro	Xilinx Virtex-II XC2V6000-ff1152-4	Xilinx Virtex-V XC5VFX70T and
			2VP50ff1517-7	(rule-oriented) and Virtex-II RC2000	Virtex-VII XC7VX485T
				(region-oriented)
	Host	processing	Not given.	1.73GHz Intel Pentium M processor	Intel Core i5 -5220 at 3 GHz, with 8
	platform			with 2GB of memory	GB of RAM.
	HDL		VHDL	Handle C	VHDL
	Experiment

Table 5 .

 5 The qualitative attributes of FPGA implementations

	Item	Biljana Petreska	Van Nguyen	Juan Quiros

Table 6 .

 6 The summing up of CUDA-GPU implementations I. The target CUDA-GPU is NVIDIA C1060 (240 cores) and the host platform is a computer with two Intel i5 Nehalem processors (8 cores).

	Item	PCUDA	PCUDASAT	TSPCUDASAT
	P system class	Cell-like (recognizer P system with	Cell-like P system with active mem-	Tissue-like P system (resolve SAT in
		active membranes)	branes (resolve SAT in linear time)	linear time)
	P system ingredi-	Elementary membranes are repre-	Each elementary is assigned a thread	Every Cell is assigned a thread block.
	ents representation	sented by thread blocks. Threads are	block. Only the objects appearing in	More objects are requested by every
		assigned to objects.	the input multisets will be assigned a	cell than the PCUDASAT simulator.
			thread.
	Operation Process	A transition step is performed in se-	1. Generation 2. Synchronization 3.	1. Generation 2. Exchange 3. Synchro-
		lection stage and execution stage.	Check-out 4. Output.	nization 4. Check-in 5. Output.
	Experiment models (a) An illustrator model to stress the	A sequential and two CUDA-based	A CUDA-based simulator.
		simulator. (b) A model aims at solv-	parallel simulators (one of them is de-
		ing a SAT problem.	signed by a hybrid method).
	Size of the simu-	512 objects and 1024 membranes for	256 objects and 4096 membranes.	256 objects and 4096 membranes.
	lated P system	(a) and 914 objects and 4096 mem-
		branes for (b).
	Speedup	7 times for (a) and 1.67 times for (b). The CUDA simulator is 63 times
			faster than the sequential simulator
			and the hybrid CUDA simulator is 9.3
			times faster than the normal one.

Table 7 .

 7 The summing up of CUDA-GPU implementations II

	Item	ABCD-GPU	ENPS-GPU	snpgpu and CuSNP
	P system class	Population Dynamics P system	Enzymatic numerical P system (cell-	Spiking Neural P system
			like)
	P system ingredi-	Environments and simulations are	Each production function and each	Spiking transition matrices, spiking
	ents representation	distributed through thread blocks.	repartition protocol element are asso-	vectors and configuration vectors are
		Rule blocks are assigned among	ciated with a thread respectively.	employed to model the target sys-
		threads.		tems.
	Operation Process	(DCBA algorithm) Selection stage (1.	1. Select the applicable programs 2.	Generate the matrix representation
		Distribution 2. Maximality 3. Proba-	Compute production functions 3. Dis-	for the model, and iterate the com-
		bility) and Execution stage.	tribute the results of the precede step	putation of next configuration until
			according to repartition protocols.	a stopping criterion is met.
	Experiment models A set of randomly generated PDP	A dummy model and an exponential	(a) SNP systems covering each com-
		systems (without biological mean-	function approximation model.	putation path sequentially, without
		ing), the ecosystem model of the		delays. (b) Extensions enabling de-
		Bearded Vulture and a tritrophic ex-		lays, support of more types of regu-
		ample model with modules defined		lar expressions and input of P-Lingua
		as features.		files.
	Size of the simu-	Running 50 simulations, 20 environ-	2 membranes, range from 1 to 120000	(a) From 1 to 16 neurons. (b) General-
	lated P system	ments and more than 20000 rule	programs.	ized sorting networks with up to 512
		blocks.		neurons.
	Speedup	7 times (Tesla C1060) compared with	6.5 times for dummy model and 10	(a) 2.31 times for the Python ver-
		the sequential model, 3 times con-	times for exponential model. 11 times	sion simulator. (b) 50 times on CUDA-
		trasted to the 4 cores CPU. The	on GeForce GTX 460M.	based version (GTX750 GPU).
		speedup for the real ecosystem, 4.9x
		with a C1060, and 18.1x using a
		Tesla K40 GPU (2880 cores). Us-
		ing the adaptative simulator with
		the tritrophic model using modules
		achieved an extra 2.5x on a K40 GPU
		and 1.8x on a P100 GPU (3584 cores).

Table 8 .

 8 The comparison of the quantitative attributes of the implementations between FPGA and CUDA-GPU

	Item	FPGA	CUDA-GPU
	Inchoate year	2003	2009
	Hardware Price	50-4500 $	150-2000 $.
	Concerned institute	Southwest Jiaotong University, Université Paris-	Universidad de Sevilla, University of the Philip-
		Est Créteil Val de Marne, Xihua University	pines Diliman
	Number of researchers	6	8
	Implemented P system type	Cell-like P system	P systems with active membranes, tissue P sys-
			tems with cell division for a SAT solution, Popula-
			tion dynamics P system, spiking Neural P system,
			enzymatic Numerical P system, evolution com-
			munication P system.
	Size of the simulated P systems	No limit for objects and membranes, but less than	No limit as long as GPU resources (mainly mem-
		1000 rules	ory) are available.
	Approximate time needed to obtain	2 months	1 month.
	a hardware system
	Speed-up effect	2×10 7 computational steps per second at present. Depends on model, implementation or GPU de-
			vice. Reported speedups from 2 to 90x compared
			to a sequential counterpart.
	Handling of non-determinism

DND algorithm and formal power serial theory. DCBA, or ad-hoc strategies (e.g. assuming confluent P systems).

Table 9 .

 9 The comparison of qualitative attributes of the implementations between FPGA and CUDA-GPUThe programming framework is flexible enough to support other ingredients, but it might require a significant revision to the parallel design.ScalabilityWith the size of the subject P systems growing, the computation clock rates decrease not too much and the hardware usages are approximately proportional to the size.Newer GPUs showed to reach higher performance, but the simulations have been demonstrated to be memory bandwidth bounded. It requires to run large models to see a positive speedup. Experiment subject Cell-like P systems with 10-1000 rules, 1 to 55 regions, 3-165 objects. (a) P systems with active membranes with 914 objects/rules and 4096 membranes, (b) SAT solutions in active membranes and cell division 256 objects/rules and 4096 membranes, (c) PDP systems up to 200000 rules on 50 simulations and 20 environments, (d) SN P systems with 16 to 512 neurons, (e) ENP systems with 2 membranes and up to 120000 programs. Difficulties

	Item	FPGA	CUDA-GPU
	Complexity of the parallel	Both region-level and system-level parallelism are re-	Both region-level and system-level parallelism are re-
	framework	alized, as well as the non-determinism (the equiprob-	alized, several strategies are carried out to sort out
		ability of the rules are guaranteed in Juan's contribu-	non-determinism (for PDP systems, based on DND al-
		tion)	gorithm in which the equiprobability is not proven
			mathematically.).
	Extensibility	Although the hardware consumptions are no more
		than 30% when rules are no more than 300, which
		means there are amount of space to accommodate new
		features of the P system, extending the existing type
		of P system to additional types needs significant revi-
		sions.

ACM Comput. Surv., Vol. 1, No. 1, Article . Publication date: June 2020.

times higher than the sequential simulator. ACM Comput. Surv., Vol. 1, No. 1, Article . Publication date: June 2020.

ACKNOWLEDGMENTS

The work of G.Z., Z. S. and S. V. is supported by the National Natural Science Foundation of China (61972324, 61672437, 61702428), Beijing Advanced Innovation Center for Intelligent Robots and Systems (2019IRS14), Artificial Intelligence Key Laboratory of Sichuan Province (2019RYJ06), the New Generation Artificial Intelligence Science and Technology Major Project of Sichuan Province (2018GZDZX0043) and the Sichuan Science and Technology Program (2018GZ0086). M.A.M-d-A., L.V-C. and M.J.P-J. also acknowledge the support of the research project TIN2017-89842-P (MABICAP), co-financed by Ministerio de Economía, Industria y Competitividad (MINECO) of Spain, through the Agencia Estatal de Investigación (AEI), and by Fondo Europeo de Desarrollo Regional (FEDER) of the European Union.

The authors are grateful to the Editor-in-Chief, Prof. Sartaj Sahni, the anonymous handling editor and all reviewers for their insightful and detailed comments on this manuscript, and are also indebted to Academician Gheorghe Păun for his useful discussions and valuable suggestions.