An Overview of Hardware Implementation of Membrane Computing Models

An Overview of Hardware Implementation of Membrane Computing Models GEXIANG ZHANG, Chengdu University of Technology, China and Southwest Jiaotong University, China ZEYI SHANG, Southwest Jiaotong University, China and Université Paris-Est Créteil Val de Marne, France SERGEY VERLAN, Université Paris-Est Créteil Val de Marne, France MIGUEL Á. MARTÍNEZ-DEL-AMOR, Universidad de Sevilla, Spain CHENGXUN YUAN, Southwest Jiaotong University, China LUIS VALENCIA-CABRERA and MARIO J. PÉREZ-JIMÉNEZ, Universidad de Sevilla, Spain 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, such as field-programmable gate arrays (FPGAs) or compute unified device architecture (CUDA)-enabled graphic processing units (GPUs). This article 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. 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.Md-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 byFondo Europeo de Desarrollo Regional (FEDER) of the European Union. Authors’ addresses: G. Zhang (corresponding author), Chengdu University of Technology, No.1, Dongsan Road, Erxianqiao, Chenghua District, Chengdu, 610059, China, Southwest Jiaotong University, 999 Xi’an Rd, Chengdu, 611756, China; email: [email protected]; Z. Shang, Southwest Jiaotong University, 999 Xi’an Rd, Chengdu, 611756, China, Université Paris-Est Créteil Val de Marne, 61 av. du général de Gaulle, Créteil, 94010, France; email: [email protected]; S. Verlan, Université Paris-Est Créteil Val de Marne, 61 av. du général de Gaulle, Créteil, 94010, France; email: [email protected]; M. Á. Martínezdel-Amor, L. Valencia-Cabrera, and M. J. Pérez-Jiménez, Universidad de Sevilla, Avda. Reina Mercedes S/N, Seville, 41012, Spain; emails: {mdelamor, lvalencia, marper}@us.es; C. Yuan, Southwest Jiaotong University, 999 Xián Rd, Chengdu, 611756, China; email: [email protected]. 1 INTRODUCTION With the arising of protocells 3.84B years ago in hydrothermal vent precipitates [26], the evolution of unicellular organisms led to the emergence of multicellularity [83]. 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 biochemical reactions, Gheorghe Păun initiated the area of membrane computing in 1998 as a theoretical com-puter science model borrowing many concepts from the cell biology [84]. The model, also called P systems, features a nested membrane-like structure delimiting regions that host objects (mod-eling cell chemicals) that 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 depend-ing on the manipulated objects, used types of rules, and the parallel execution strategy (called derivation mode). 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 systems, (5) dealing with multisets, (6) evolving by rewriting-like rules, then P systems are obligatory/unavoidable and the unique models dealing with all these features [85]. 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 [76]. 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 Reference [87] summarizing 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 [25, 114], meta-heuristic algorithms for optimization problems [116, 117, 119], robot controllers [14, 109], path planning [79, 81, 110], and fault diagnosis for electrical power systems [90, 107, 108]. Due to its historical background, membrane computing was also used as a modeling framework for biological and ecological subjects such as artificial life [94], photosynthesis [75], p53 protein signaling pathway [95], myxobacterial colony [68], and biopolymer duplication [56]. (Please see a recent overview in Reference [118]). For 20 years many applications of P systems have arisen, not only in Theoretical Computer Science, but also in Computational Modelling, Robotics, Optimization, and so on, 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 [112]. Therefore, there is a need for simulators to build model validation tools, assistants for formal verification, environments for virtual experimentation, model calibration tools, and so on. So far, most of the simulation software aims at reproducing only the computation results instead of procedures. However, there has also been 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 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 [12]. Natural Computing models can provide solutions where the instructions are close 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, analyzing all the challenges and current solutions. In addition, a number of these applications mentioned above only theoretically profit from the potential speedup promised by the model, as there are no truly parallel implementations available. However, 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) [51,54]. As pointed out in Reference [43] “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 transistors in chips cannot follow this law, at least it cannot be doubled within two years [2]. With transistors shrunk to nanoscale, quantum effects stand out [66] 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 [111]. 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 [1]. Parallel computing has the potential to further uplift computing power provided that the density of transistor is constant [101] with multicore and multithread architecture, although heat dissipation and interconnect issues would be challenges [8]. 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 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 [118]. 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 such as creation and dissolution of inner membranes, the self-replication or autopoiesis [29] of the whole cell-like entirety that works as a computing unit, the communication by symport and antiport of objects [4], and so on. 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 [85]. 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 [22]. 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, and so on. 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 [88]. 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 NP-complete [21]. Known algorithms and heuristics do not parallelize well, so special heuristics were developed to quickly compute the desired multiset of rules. We decided to present these heuristics uniformly in terms of multi-criteria optimization (see Sections 4and 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 article is organized as follows: Section 2gives a brief description of the capabilities and of the architecture of hardware used, Section 3recalls the definition of the most general model of P systems, based on the formal framework [35]. Next, Section 4expresses 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,and8give more details on existing simulation approaches. At the end, conclusions and future research directions are discussed. 2 THE SELECTION OF HARDWARE 2.1 FPGA Hardware 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. FPGA is developed by 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 [67]). From structural point of view, an FPGA is an array of configurable logic blocks (CLB) inlaid in the matrix of interconnects. CLBs comprise slice-organized logic cells which are arranged in a way named \emph{look-up table} (LUT). LUTs are used to implement different combinatorial circuits such as basic gates, decoders, encoders, and multiplexers. Logic cells contain a type of storage element called flip-flop (FF) which is used as memory elements in sequential circuits. To perform particular operations, 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. 2.2 CUDA-enabled GPU Hardware Nowadays, multi-core architecture CPU is the mainstream. However, 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) [52]. 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 be concerned 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 that can handle other computing assignments except for graphics processing. This predicament has changed for the arise of compute unified device architecture, known as CUDA, from 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 a universal parallel computing device that 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. The CUDA-enabled GPU does not work alone, but form a heterogeneous computing architecture with the CPU, 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 [52]. 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 [52]. Although current GPUs contain cache memories, 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 the same block. Threads also have fast access to their own registers 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) that 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 that an SM fulfills in its executions. An SM can manage multiple warps that 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 [61]. 2.3 Other Hardware There were attempts to simulate certain types of P systems on micro-processor–based architectures [47]. 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 [93] 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 article. 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 [28, 103]. 3 THE MODEL OF P SYSTEMS We suppose the reader has a knowledge of basic notions from formal language theory and membrane computing. We refer to References [37, 87, 91] for missing details. We will also closely follow References [35, 105] 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 Mwill be represented by a string where the number of occurrences of each letter corresponds to its multiplicity in M. We will denote by |M|thesizeofthemultisetMand by |M|athe number of elements ain M. There exist many variants of P systems (see, e.g., Reference [85]). In this article the presentation will be given in terms of the formal framework for P systems introduced in Reference [35], see also References [104,105]. 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 four 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 Reference [34], but corresponding definitions are too complex for the purpose of this article. 3.1 Network of Cells As pointed out in References [33,35,104,105], 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 Reference [35]. 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. References [35,104,105]definesome 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 Reference [35]. We remark that the definition from Reference [34] is slightly different, however, both models coincide when the structure of the system does not evolve. Definition 3.1 ([35]). Anetwork of cells of degree n≥1 is a construct Π=(n,V,w,Inf ,R), where (1) nis the number of cells; (2) Vis an alphabet; (3) w=(w1,...,wn)where wi∈V◦, for all 1 ≤i≤n,isthefinite multiset initially associated to cell i; (4) Inf =(Inf1,...,Infn),whereInfi⊆V, for all 1 ≤i≤n,istheset of symbols occurring infinitely often in cell i(in most of the cases, only one cell, called the environment, will contain symbols occurring with infinite multiplicity); (5) Ris a finite set of rules of the form (X→Y;P,Q), where X=(x1,...,xn),Y=(y1,...,yn),xi,yi∈V◦,1≤i≤n, are vectors of multisets over Vand P=(p1,...,pn),Q=(q1,...,qn),pi,qi,1≤i≤nare finite sets of multisets over V. We will also use the notation (omitting pi,qi,xior yiif they are empty) (1,x1)...(n,xn)→(1,y1)...(n,yn);[(1,p1)...(1,pn)]; [(1,q1)...(n,qn)]. The above rule is applied as follows: objects xifrom cellsiare rewritten into objectsyjproduced in cells j,1≤i,j≤n, if every cell k,1≤k≤n, contains all multisets from pkand does not contain any multiset from qk. 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 Cof Πis defined as an n-tuple of multisets over V(u1,...,un)satisfying ui∩Infi=∅,1≤i≤n. To define the computation in network of cells according to some derivation modeδthe following functions should be specified: •Applicable (Π,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. •Apply (Π,C,R)– the function allowing to compute the configuration obtained by the parallel application of the multiset of rules Rto the configuration C. •Halt(Π,C,δ)– a predicate that yields true if Cis a halting configuration of the system Π (in some derivation mode δ). •Result(Π,C)– a function giving the result of the computation of the P system Πwhen the halting configuration Chas been reached. Then the computation is a sequence of transitions where each transition step C⇒Cis defined as C=Apply (Π,C,R),for some R∈Applicable (Π,C,δ). As usual, this sequence starts with the initial configuration and ends with the final configuration for which the halting predicate Halt yields true. In more formal terms, the result of the computation of a network of cells Π=(n,V,w,Inf ,R)working in the derivation mode δis defined as follows (we refer to Reference [35] for more technical details): Result(Π) ={Result(Π,z):w⇒∗z,Halt(Π,z,δ)=true and Halt(Π,x,δ)=falsefor any x:w⇒∗x⇒+z}. We remark that Reference [35] gives an algorithm to compute the set of multisets of applicable rules, denoted as Applicable (Π,C, asyn), the algorithm to compute Apply(Π, C, R), and several definitions for Halt and R esult 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 Applicable (Π,C, asyn) 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 (max), which is defined as follows: Applicable (Π,C,max ) = {R ⊆ Applicable (Π,C, asyn) | R ∈ Applicable (Π,C, asyn) : R  R}. 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 Reference [105] for a description of several other derivation modes. Example 3.2. Consider the system Π=(O,w1,R), having the alphabet O={a,b,c}and the set of rules R={r1:(1,ab)→(1,abc);r2:(1,bbc)→(1,abb)}. Consider the configuration C= (1,a3b4c2).Then,Applicable (Π,C,asyn)=r1,r2 1,r3 1,r2,r2 2,r1r2,r2 1r2. The maximally parallel set of multisets of rules is the following: Applicable (Π,C,max)=r3 1,r2 2,r2 1r2. The application of the multiset of rules r2 1r2on Cyields (1,a4b4c3). Example 3.3. Consider the system Π=(O,w1,w2,R), with O={a,b,c}and R={r1: (1,a)(2,b)→(1,b)(2,a);r2:(1,b)→(2,b)}. Consider the configuration C=(1,a4)(2,b).It is easy to observe that at each step only a single rule is applicable (r1or r2alternatively). The only possible sequence of rule applications is (r1r2)4after which no rule is applicable anymore. This sequence moves all symbols afrom cell 1 to cell 2. We remark that the above rules do not rewrite objects, but only move them in the structure. 3.2 Examples of P Systems Here, we give the informal description of some main variants of P systems that were targeted for an implementation. Transitional (cell-like) P systems [77]. This model considers that cells are organized in a tree structure and uses rules of following type: (i,u)→(i,u)(j,u)(k1,v1)...(km,vm),wherejis the parent of iand iis the parent of k1,...,km,m≥0. When |u|>1 corresponding rules are called cooperative. Symport/antiport P systems [4]. 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 (i,u)(j,v)→ (i,v)(j,u)or (i,u)→(j,u). Populational Dynamics P (PDP) systems [13,92]. 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: (i,u)(j,v)→(i,u)(j,v),whereiis parent of jand (2) communication between tree roots: (ri,x)→(r1,x1)...(rk,xk),wherer1,...,rkare 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 [24,49]. 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 (i,ak)→(k1,ax1)...(km,axm);(i,E),whereE is the regular expression checking for the contents of cell i,ais the single symbol used in the alphabet and k1,...,kmare the cells linked to cell i. 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 [78,86]. A multiset Mcan also be seen as a function M:V→ N. The model of numerical P systems extends the notion of multiset to M:V→Rn. Hence, the 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 References [41, 44, 57, 60], hardware architectures aiming at parallel processing and communication, and the application of rules are developed. In Reference [11], a formal exposition of nondeterministic 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 Applicable (Π,C, δ ) (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 Reference [71] 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 [31]. Before the first correct solution is obtained, some invalid solutions should be rejected. This approach is called indirect straightforward approach [71]. 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 References [71, 73] 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 References [30, 39, 46, 48, 98, 99]. 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 Reference [71], 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 large number of cases. As shown in References [88, 106], if at each step the number of solutions can be expressed as 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 Reference [106] 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 Reference [71]. A similar algorithm can also be found in References [38,40]. This algorithm works in two 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 and their applicability is increased up to the maximal still possible value, except that the last rule in the first phase keeps its original applicability value. A variant of DND, named DND-P, became popular in the simulation of Population Dynamics P (PDP) systems [65]. Together with another algorithm, Direct distribution based on Consistent Blocks Algorithm (DCBA) [63], it was employed for the engine of the PDP system simulator on CUDA [62]. Non-determinism. One of the difficulties of the above approaches is the handling of 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 (Applicable (Π,C,δ)). 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 that 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 References [88, 106]. However, the corresponding implementation is limited to some particular derivation modes and particular classes of P systems. 6 AN OVERVIEW OF EXISTING FPGA SIMULATIONS With the advent of reconfigurable hardware that 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 [82]. Since then, two simulation approaches emerged, considering regions or rules as basic processing units. 6.1 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) [82]. They could realize several interesting features such as 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, 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 u → v (v1, ini )(v2, out ),wherev1 is the string to be sent into lower-immediate membrane labeled i, v2 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 u ≤ w, this rule is applicable and this module will generate a signal Applicable = 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, v1 and v2 are stored in corresponding registers. 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 terminate 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. 6.1.2 Nguyen Simulation. In this implementation, a parallel computing platform simulating membrane computing based on FPGA named Reconfig-P is developed [69,74]. 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 [69]. 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. 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 immediately after an applicable is applied to consume some objects. 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 Reference [69] two solution strategies, the space-oriented strategy and the time-oriented strategy, are proposed. Tables 1and 2summarize the different strategies of two resource conflict resolutions and their modifications in the rule-based and region-based design [70,72]. 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 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. 6.2 Rule-based Simulations Rule-based approaches consider evolution rules as processing units performing the update of multiplicities and of membrane structures. 6.2.1 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. In the design phase, a processing unit corresponds to a potential infinite while loop that contains the Handel C codes related to the rule application. The information associated to execution and Table 1. The Comparison of the Time-oriented and Space-oriented Conflict Resolution Strategy Resource conflict resolution Time-oriented strategy (a) Construct a conflict matrix in which each row is a quadruple (p,q,r,s). pis the object competed by multiple rules. qis the region where pis produced or consumed. ris the set of the conflicting rules, sis the size of set r. (b) Insert delay statement among conflict rules such that the updating operations of multiplicity of pcan be executed in distinct clock cycles. The number of delays is equal to s. 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 s 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 run-time. The role of the multiset replication coordinator in rule-oriented 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 that send objects to the considered region to write the new values of the objects competed by multiple rules. 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. If there are objects traversing regions that 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. 6.2.2 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 to not compute the complete solution, i.e., Applicable (Π,C, δ ) 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 [88, 106]. Devising an algorithm that 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 Table 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 containment relationships 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 ruleoriented 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 1bit 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. 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: NBVariants(Π,C,δ), which gives the cardinality of Applicable (Π,C,δ),and Variant(n,Π,C,δ), which returns the nth element of Applicable (Π,C,δ). 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 r1:ab →uand r2:bc →v. Assume that the derivation mode is maximal parallelism (max). Suppose that Na,Nb,and Ncrepresent the number of objects a,b,and cin C.LetN1=min(Na,Nb),N2=min(Nb,Nc), N=min(N1,N2),ki=NiN,1≤i≤2, where denotes the positive subtraction. Let also p,q=0,1,2,...,N. From the dependency graph, we can deduce the following: Applicable (Π,C,max)= p+q=N rp+k1 1rq+k2 2 NBVariants(Π,C,max)=N+1 Variant(n,Π,C,max)=rN−n+1+k1 1rn−1+k2 2. Example 6.1. Consider a configuration where Na=5, Nb=5,and Nc=3. It can be easily verified that N1=min(5,5)=5, N2=min(5,3)=3, N=min(5,3)=3, k1=N1N=5−3=2, k2=N2N=3−3=0. Hence, we can enumerate the elements of Applicable (Π,C,max)as below: Applicable (Π,C,max)1={r3+2 1r0+0 2,r2+2 1r1+0 2,r1+2 1r2+0 2,r0+2 1r3+0 2}={r5 1,r4 1r2,r3 1r2 2,r2 1r3 2}. The same result can be easily obtained by using formal power series associated to contextfree languages. In this case, any maximal rule combination is a part of the language LN = {r1 pr2 q | p + q = N }. It is quite easy to observe that the number of words of length N in LN is exactly the same as the number of words of same length in the language L = {r1∗r2∗}. This last language is regular and its generating function is q0 (x ) = 1/(1 − x )2.Thenth coefficient of the expansion of q0 (x ) is equal to n + 1([xn ]q0 = n + 1), which immediately gives NBVariants(Π,C,max ) = N + 1. The Variant (n, Π, c,max ) is computed using an algorithm that performs a weighted breadth-first search of the decomposition of n with respect to the number of variants found on each branch of the execution of automaton for L. 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 Reference [106] exhibited a speedup of order 105. 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 NBVariants(Π,C,max ), Variant (n, Π,C,max ) 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. An independent stage computes the maximal instance of each rule by means of the dividing operation and MIN logic operation. Assignment stage is 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. In order to more clearly show the FPGA implementation methods of P system proposed by the above three research groups, their methods are summarized and compared from the quantitative and qualitative perspectives in Table 4 and Table 5. Table 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 Lausanne (EPFL) University of South Australia University of Seville Target FPGA Xilinx Virtex-II Pro 2VP50ff1517-7 Xilinx Virtex-II XC2V6000ff1152-4 (rule-oriented) and Virtex-II RC2000 (regionoriented) Xilinx Virtex-V XC5VFX70T and Virtex-VII XC7VX485T Host processing platform Not given. 1.73 GHz Intel Pentium M processor with 2 GB of memory Intel Core i5—5220 at 3 GHz, with 8 GB of RAM. HDL VHDL Handle C VHDL Experiment Subjects Cell-like P systems with following characteristics: membrane dissolution and creation, objects exchange between upperand 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 region-oriented design, the objects are cell-like P systems containing hierarchical regions and tissue-like connected regions. The rule range is [10,50], [1,25] for regions, [3,200] for objects. The extent of the inter-region communication is [80,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 four subjects differ in the rule dependencies that form chains: circular, 2-circular, linear, opposite. The object number range is [10,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×105 to 1000 ×105. The hardware consumption extent is 1.55% to 21.43% (LUTs). For regionoriented 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 that can be executed is a P system with 550 rules, 1,280 communication channels, 1,100 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 ×107computational steps per second. The run-time for each experiment subject ranges from 3.017 × 10−5sto 4.174 ×10−5s. 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. 7 AN OVERVIEW OF EXISTING CUDA SIMULATIONS In Reference [61], 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; Table 5. The Qualitative Attributes of FPGA Implementations Item Biljana Petreska Van Nguyen Juan Quiros 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 is 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 right-hand 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 regionoriented design, they are implicitly expounded through integrating them into the region processing units. The logic of rules is distributed along 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 rises 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 nondeterminism. Drawback Partially parallelism, no nondeterminism The equiprobability of DND algorithm has not been proven theoretically. The types of P systems that can be implemented are confined to those whose rules can be represented as regular language. •A medium scalability degree: The amount of resources depend on the GPU model, e.g., a K40 includes 2,880 cores and 12 GBytes of memory. If the resources of a GPU are 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, and generating the right-hand side. This strategy is necessary 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. 7.1 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 [18]. This simulator performs only one computation out of the whole tree 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 constraints, 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 [18], 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 requires 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 [18]. This way, the simulator can choose any path in the computation tree, since the aim is to find a halting configuration. The GPU simulator takes 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, 4 GB memory) by using two benchmarks [61]: a simple test P system designed in a convenient manner (up to 7× of speedup), and a family of P systems designed to solve different instances of the SAT problem (1.67×of speedup). 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 venture. Slower than FPGA. Sometimes big concessions are indispensable, as the architectures are unchangeable. Contributions Introducing the reconfigurable hardware to develop the real parallel architectures that 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 that implement the target models as approximate as possible on FPGA devices. The expense we should pay is also considerablly 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. ACKNOWLEDGMENTS 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. REFERENCES [1] Rick Merritt. 2017. Roadmap Says CMOS Ends ∼2024. IRDS points to chip stacks, new architectures. Retrieved from https://web.archive.org/web/20170324022546/,https://www.eetimes.com/document.asp?doc_id=1331517. [2] Gordon Moore. 1975. Progress in digital integrated electronics. International Electron Devices Meeting. IEEE, 11–13. [3] Oana Agrigoroaiei, Gabriel Ciobanu, and Andreas Resios. 2010. Evolving by maximizing the number of rules: Complexity study. In Proceedings of the 10th International Workshop on Membrane Computing (WMC’09) (LNCS),Gheorghe Păun, Mario J. Pérez-Jiménez, Agustín Riscos-Núñez, Grzegorz Rozenberg, and Arto Salomaa (Eds.). Springer, 149–157. [4] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter. 2002. Molecular Biology of the Cell (4th ed.). Garland. [5] Artiom Alhazov. 2005. Maximally parallel multiset-rewriting systems: Browsing the configurations. In Proceedings of the 3rd Brainstorming Week on Membrane Computing (RGNC Report), M. A. Gutiérrez-Naranjo, A. Riscos-Núñez, F. J. Romero-Campero, and D. Sburlan (Eds.). 1–10. [6] Santiago Alonso, Luis Fernández, Fernando Arroyo, and Javier Gil. 2008. A circuit implementing massive parallelism in transition P systems. Int. J. Inform. Technol. Knowl. 2, 1 (2008), 35–42. [7] Santiago Alonso, Luis Fernández, Fernando Arroyo, and Javier Gil. 2008. Main modules design for a HW implementation of massive parallelism in transition P-systems. Artif. Life Robot. 13, 1 (2008), 107–111. [8] Himanshu Kushwah Anil Sethi. 2015. Multicore processor technology—Advantages and challenges. Int. J. Res. Eng. Technol. 4, 9 (2015), 87–89. [9] Alberto Arteta, Luis Fernández, and Javier Gil. 2008. Algorithm for application of evolution rules based on linear diofantic equations. In Proceedings of the 10th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC’08), Viorel Negru, Tudor Jebelean, Dana Petcu, and Daniela Zaharie (Eds.). IEEE Computer Society, 496–500. [10] Jambhlekar P. Arun, Manoj Mishra, and Sheshasayee V. Subramaniam. 2011. Parallel implementation of MOPSO on GPU using OpenCL and CUDA. In Proceedings of the 18th International Conference on High Performance Computing. 1–10. [11] Angel V. Baranda, Fernando Arroyo, Juan Castellanos, and Rafael Gonzalo. 2001. Towards an electronic implementation of membrane computing: A formal description of non-deterministic evolution in transition P systems. In Proceedings of the 7th International Workshop on DNA-Based Computers: DNA Computing (LNCS), N. Jonoska and N. C. Seeman (Eds.), Vol. 2340. Springer, 350–359. [12] Sankar Basu, Randal E. Bryant, Giovanni De Micheli, Thomas Theis, and Lloyd Whitman. 2019. iFLEX: A fully open-source high-density field-programmable gate array (FPGA)-based hardware co-processor for vector similarity searching. Proc. IEEE 107, 1 (2019), 11–18. [13] Francesco Bernardini and Marian Gheorghe. 2004. Population P systems. J. Univ. Comput. Sci. 10, 5 (2004), 509–539. [14] Catalin Buiu, Cristian Vasile, and Octavian Arsene. 2012. Development of membrane controllers for mobile robots. Inf. Sci. 187 (2012), 33–51. [15] Francis G. C. Cabarle, Henry N. Adorna, Miguel A. Martínez-del-Amor, and Mario J. Pérez-Jiménez. 2012. Improving GPU simulations of spiking neural P systems. Roman. J. Inf. Sci. Technol. 15, 1 (2012), 5–20. [16] Jym P. Carandang, Francis G. C. Cabarle, Henry N. Adorna, Nestine H. S. Hernandez, and Miguel Á. Martínez-delAmor. 2019. Handling non-determinism in spiking neural P systems: Algorithms and simulations. Fundam. Inform. 164 (2019), 139–155. [17] Jym P. A. Carandang, John M. B. Villaflores, Francis G. C. Cabarle, Henry N. Adorna, and Miguel A. Martínez-delAmor. 2017. CuSNP: Spiking neural P systems simulators in CUDA. Roman. J. Inform. Sci. Technol. 20, 1 (2017), 57–70. [18] José M. Cecilia, José M. García, Ginés D. Guerrero, Miguel A. Martínez-del-Amor, Ignacio Pérez-Hurtado, and Mario J. Pérez-Jiménez. 2010. Simulation of P systems with active membranes on CUDA. Brief. Bioinform. 11, 3 (2010), 313–322. [19] José M. Cecilia, José M. García, Ginés D. Guerrero, Miguel A. Martínez-del-Amor, Mario J. Pérez-Jiménez, and Manuel Ujaldon. 2012. The GPU on the simulation of cellular computing models. Soft Comput. 16, 2 (2012), 231–246. [20] Gabriel Ciobanu, Solomon Marcus, and Gheorghe Păun. 2009. New strategies of using the rules of a P system in a maximal way: Power and complexity. Roman. J. Inform. Sci. Technol. 12, 2 (2009), 157–173. [21] Gabriel Ciobanu and Andreas Resios. 2009. Complexity of evolution in maximum cooperative P systems. Nat. Comput. 8, 4 (31 Jan. 2009), 807. [22] Gabriel Ciobanu and Guo Wenyuan. 2003. P systems running on a cluster of computers. In Proceedings of the International Workshop Membrane Computing (WMC’03) (Lecture Notes in Computer Science), Carlos Martín-Vide, Giancarlo Mauri, Gheorghe Páun, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 2933. Springer, 123–139. [23] Daniele D’Agostino, Giulia Pasquale, and Ivan Merelli. 2014. A fine-grained CUDA implementation of the multiobjective evolutionary approach NSGA-II: Potential impact for computational and systems biology applications. In Proceedings of the 11th International Meeting of Computational Intelligence Methods for Bioinformatics and Biostatistics (CIBB’14) (LNCS), Clelia Di Serio, Pietro Liò, Alessandro Nonis, and Roberto Tagliaferri (Eds.), Vol. 8623. Springer, 273–284. [24] Ren T. A. de la Cruz, Francis G. C. Cabarle, and Henry N. Adorna. 2019. Generating context-free languages using spiking neural P systems with structural plasticity. J. Memb. Comput. 1, 3 (2019), 161–177. [25] Daniel Díaz-Pernil, Miguel A. Gutiérrez-Naranjo, and Hong Peng. 2019. Membrane computing and image processing: A short survey. J. Memb. Comput. (04 Feb. 2019). [26] Matthew S. Dodd, Dominic Papineau, Tor Grenne, John F. Slack, Martin Rittner, Franco Pirajno, Jonathan O’Neil, and Crispin T. S. Little. 2017. Evidence for early life in Earth’s oldest hydrothermal vent precipitates. Nature 543, 7643 (2017), 60–64. [27] Matthias Ehrgott. 2005. Multicriteria Optimization (2nd ed.). Springer. [28] H. Esmaeilzadeh, A. Sampson, L. Ceze, and D. Burger. 2012. Neural acceleration for general-purpose approximate programs. In Proceedings of the 45th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO’12). IEEE Computer Society, 449–460. [29] R. Uribe, F. Varela, and H. Maturana. 1974. Autopoiesis: The organization of living systems, its characterization and a model. BioSystems 5, 4 (1974), 187–196. [30] Luis Fernández, Fernando Arroyo, Ivan Garcia, and Gines Bravo. 2007. Decision trees for applicability of evolution rules in transition P systems. Inf. Theor. Applic. 14, 3 (2007), 223–230. [31] Luis Fernández, Fernando Arroyo, Jorge A. Tejedo, and Juan Castellanos. 2006. Massively parallel algorithm for evolution rules application in transition P systems. In Proceedings of the 7th Workshop on Membrane Computing, Hendrik Jan Hoogeboom, Gheorghe Păun, and Grzegorz Rozenberg (Eds.). Universiteit Leiden, 337–343. [32] Luis Fernández, Victor J. Martínez, Fernando Arroyo, and Luis F. Mingo. 2005. A hardware circuit for selecting active rules in transition P systems. In Proceedings of the 7th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC’05), Daniela Zaharie, Dana Petcu, Viorel Negru, Tudor Jebelean, Gabriel Ciobanu, Alexandru Cicortas, Ajith Abraham, and Marcin Paprzycki (Eds.). IEEE Computer Society, 415–418. [33] Rudolf Freund, Alberto Leporati, Giancarlo Mauri, Antonio E. Porreca, Sergey Verlan, and Claudio Zandron. 2013. Flattening in (Tissue) P systems. In Proceedings of the 14th International Conference on Membrane Computing (CMC’13) (Lecture Notes in Computer Science), Artiom Alhazov, Svetlana Cojocaru, Marian Gheorghe, Yurii Rogozhin, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 8340. Springer, 173–188. [34] Rudolf Freund, Ignacio Pérez-Hurtado, Agustín Riscos-Núñez, and Sergey Verlan. 2013. A formalization of membrane systems with dynamically evolving structures. Int. J. Comput. Math. 90, 4 (2013), 801–815. [35] Rudolf Freund and Sergey Verlan. 2007. A formal framework for static (Tissue) P systems. In Proceedings of the 8th International Workshop on Membrane Computing (WMC’07) (Lecture Notes in Computer Science), George Eleftherakis, Petros Kefalas, Gheorghe Paun, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 4860. Springer, 271–284. [36] Zsolt Gazdag and Gábor Kolonits. 2019. A new method to simulate restricted variants of polarizationless P systems with active membranes. J. Memb. Comput. 1, 4 (2019), 251–261. [37] Marian Gheorghe, Andrei Păun, Sergey Verlan, and Gexiang Zhang. 2017. Membrane Computing, Power and Complexity. Springer Berlin, 1–16. [38] Francisco J. Gil, Luis Fernández, Fernando Arroyo, and Juan Alberto de Frutos. 2008. Parallel algorithm for P systems implementation in multiprocessors. In Proceedings of the 13th International Symposium on Artificial Life and Robotics (AROB’08), M. Sugisaka and H. Tanaka (Eds.). 10–25. [39] Francisco J. Gil, Luis Fernández, Fernando Arroyo, and Jorge A. Tejedor. 2007. Delimited massively parallel algorithm based on rules elimination for application of active rules in transition P systems. In Proceedings of the 5th International Conference on Information Research and Applications (i.TECH’07), Krassimir Markov and Krassimira Ivanova (Eds.), Vol. 1. Institute of Information Theories and Applications FOI ITHEA, Bulgaria, 182–188. [40] Francisco J. Gil, Jorge A. Tejedor, and Luis Fernández. 2008. Fast linear algorithm for active rules application in transition P systems. In Algorithmic and Mathematical Foundations of the Artificial Intelligence (International Book Series INFORMATION SCIENCE & COMPUTING), Krassimir Markov, Krassimira Ivanova, and Ilia Mitov (Eds.), Vol. Supple. Institute of Information Theories and Applications FOI ITHEA, Sofia, Bulgaria, 35–44. [41] Fernando Arroyo Ginés Bravo, Luis Fernández, and Juan Frutos. 2008. A hierarchical architecture with parallel communication for implementing P systems. Inf. Technol. Knowl. 2, 1 (2008), 43–48. [42] Sandra M. Gomez-Canaval, Abraham Gutiérrez, and Santiago Alonso. 2008. Hardware implementation of P systems using microcontrollers. An operating environment for implementing a partially parallel distributed architecture. In Proceedings of the 10th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC’08), Viorel Negru, Tudor Jebelean, and Dana Petcuand Daniela Zaharie (Eds.). IEEE, 489–495. [43] Ananth Grama, George Karypis, Vipin Kumar, and Anshul Gupta. 2003. Introduction to Parallel Computing (2nd ed.). Addison-Wesley. [44] Abraham Gutiérrez, Luis Fernández, Fernando Arroyo, and Santiago Alonso. 2007. Hardware and software architecture for implementing membrane systems: A case of study to transition P systems. In Proceedings of the 13th International Meeting on DNA Computing (LNCS), Max H. Garzon and Hao Yan (Eds.), Vol. 4848. Springer, 211–220. [45] Abraham Gutiérrez, Luis Fernández, Fernando Arroyo, and Santiago Alonso. 2008. Suitability of using microcontrollers in implementing new P-system communications architectures. Artif. Life Robot. 13, 1 (01 Dec. 2008), 102–106. [46] Abraham Gutiérrez, Luis Fernández, Fernando Arroyo, and Ginés Bravo. 2007. Optimizing membrane system implementation with multisets and evolution rules compression. In Proceedings of the 8th Workshop on Membrane Computing, G. Ekeftherakis, P. Kefalas, and Gh. Păun (Eds.). 345–362. [47] Abraham Gutiérrez, Luis Fernández, Fernando Arroyo, and Victor J. Martínez. 2006. Design of a hardware architecture based on microcontrollers for the implementation of membrane systems. In Proceedings of the 8th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing, Viorel Negru, Dana Petcu, Daniela Zaharie, Ajith Abraham, Bruno Buchberger, Alexandru Cicortas, Dorian Gorgan, and Joël Quinqueton (Eds.). IEEE Computer Society, 350–353. [48] Francisco Javier Gil, Luis Fernández, Fernando Arroyo, and Jorge Tejedor. 2008. Delimited massively parallel algorithm based on rules elimination for application of active rules in transition P systems. Inf. Technol. Knowl. 2, 1 (2008), 56–61. [49] Z. B. Jimenez, Francis G. C. Cabarle, Ren T. A. de la Cruz, Kelvin C. Buño, Henry N. Adorna, Nestine H. S. Hernandez, and Xiangxiang Zeng. 2019. Matrix representation and simulation algorithm of spiking neural P systems with structural plasticity. J. Memb. Comput. 1, 3 (2019), 145–160. [50] Ignacy Kaliszewski, Janusz Miroforidis, and Dmitry Podkopaev. 2016. Multiple Criteria Decision Making by Multiobjective Optimization–A Toolbox. Int. Series in Operations Research & Management Science, Vol. 242. Springer. [51] Takahiro Katagiri. 2019. High-performance computing basics. In The Art of High Performance Computing for Computational Science, Vol. 1, Techniques of Speedup and Parallelization for General Purposes, Masaaki Geshi (Ed.). Springer, 1–25. [52] David Kirk and Wen-Mei Hwu. 2010. Programming Massively Parallel Processors: A Hands On Approach.Morgan Kaufmann. [53] Julien Legriel. 2011. Multi-Criteria Optimization and Its Application to Multi-Processor Embedded Systems. Universite de Grenoble, Grenoble, France. [54] Jianping Kelvin Li, Jia-Kai Chou, and Kwan-Liu Ma. 2015. High performance heterogeneous computing for collaborative visual analysis. In Proceedings of the SIGGRAPH Asia Visualization in High Performance Computing Conference. ACM, 12:1–12:4. [55] Luis F. Macías-Ramos, Miguel A. Martínez-del-Amor, and Mario J. Pérez-Jiménez. 2015. Simulating FRSN P systems with real numbers in P-Lingua on sequential and CUDA platforms. In Proceedings of the 16th International Conference on Membrane Computing (CMC’15) (LNCS), G. Rozenberg, A. Salomaa, J. M. Sempere, and C. Zandron (Eds.). 262–276. [56] Vincenzo Manca. 2019. From biopolymer duplication to membrane duplication and beyond. J. Memb. Comput. 1, 4 (2019), 292–303. [57] Víctor Martínez, Santiago Alonso, and Abraham Gutiérrez. 2010. Hardware circuit for the application of evolution rules in a transition P-system. Artif. Life Robot. 15, 1 (01 Aug. 2010), 89–92. [58] Victor Martínez, Fernando Arroyo, Abraham Gutiérrez, and Luis Fernández. 2006. Hardware implementation of a bounded algorithm for application of rules in a transition P-system. In Proceedings of the 8th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC’06), Viorel Negru, Dana Petcu, Daniela Zaharie, Ajith Abraham, Bruno Buchberger, Alexandru Cicortas, Dorian Gorgan, and Joel Quinqueton (Eds.). IEEE, 343–349. [59] Victor Martínez, Luis Fernández, Fernando Arroyo, and Abraham Gutiérrez. 2007. HW implementation of a optimized algorithm for the application of active rules in a transition P-system. Inf. Theor. Applic. 14, 4 (2007), 324–331. [60] Victor J. Martínez, Fernando Arroyo, Abraham Gutiérrez, and Luis Fernández. 2006. Hardware implementation of a bounded algorithm for application of rules in a transition P-system. In Proceedings of the 8th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC’06), Viorel Negru, Dana Petcu, Daniela Zaharie, Ajith Abraham, Bruno Buchberger, Alexandru Cicortas, Dorian Gorgan, and Joël Quinqueton (Eds.). IEEE Computer Society, 343–349. [61] Miguel A. Martínez-del-Amor, Manuel García-Quismondo, Luis F. Macías-Ramos, Luis Valencia-Cabrera, Agustin Riscos-Núñez, and Mario J. Pérez-Jiménez. 2015. Simulating P systems on GPU devices: A survey. Fundam. Inform. 136, 3 (2015), 269–284. [62] Miguel A. Martínez-del-Amor, Luis F. Macías-Ramos, Luis Valencia-Cabrera, and Mario J. Pérez-Jiménez. 2016. Parallel simulation of population dynamics P systems: Updates and roadmap. Nat. Comput. 15, 4 (2016), 565–573. [63] Miguel A. Martínez-del-Amor, Ignacio Pérez-Hurtado, Manuel García-Quismondo, Luis F. Macías-Ramos, Luis Valencia-Cabrera, Álvaro Romero Jiménez, Carmen Graciani Díaz, Agustin Riscos-Núñez, Maria Angels Colomer, and Mario J. Pérez-Jiménez. 2012. DCBA: Simulating population dynamics P systems with proportional object distribution. In Proceedings of the 13th International Conference on Membrane Computing (CMC’12) (Lecture Notes in Computer Science), Erzsébet Csuhaj-Varjú, Marian Gheorghe, Grzegorz Rozenberg, Arto Salomaa, and György Vaszil (Eds.). 257–276. [64] Miguel Á. Martínez-del-Amor, Ignacio Pérez-Hurtado, David Orellana-Martín, and Mario J. Pérez-Jiménez. 2020. Adaptative parallel simulators for bioinspired computing models. Fut. Gen. Comput. Syst. 107 (2020), 469–484. [65] Miguel A. Martínez-del-Amor, Ignacio Pérez-Hurtado, Mario J. Pérez-Jiménez, Agustin Riscos-Núñez, and M. Angels Colomer. 2010. A new simulation algorithm for multienvironment probabilistic P systems. In Proceedings of the 5th International Conference on Bio-Inspired Computing: Theories and Applications,M.Gong,L.Pan,T.Song,andG. Zhang (Eds.). 59–68. [66] Neil Mathur. 2002. Beyond the silicon roadmap. Nature 419 (Oct. 2002), 573–575. [67] George H. Mealy. 1955. A method for synthesizing sequential circuits. Bell Syst. Tech. J. 34, 5 (1955), 1045–1079. [68] Anthony Nash and Sara Kalvala. 2019. A P system model of swarming and aggregation in a Myxobacterial colony. J. Memb. Comput. 1, 2 (2019), 103–111. [69] Van Nguyen. 2010. An Implementation of the Parallelism, Distribution and Nondeterminism of Membrane Computing Models on Reconfigurable Hardware. Ph.D. Dissertation. University of South Australia. [70] Van Nguyen, David Kearney, and Gianpaolo Gioiosa. 2007. Balancing performance, flexibility, and scalability in a parallel computing platform for membrane computing applications. In Proceedings of the 8th International Workshop on Membrane Computing (LNCS), G. Eleftherakis, P. Kefalas, Gh. Păun, G. Rozenberg, and A. Salomaa (Eds.), Vol. 4860. Springer, 385–413. [71] Van Nguyen, David Kearney, and Gianpaolo Gioiosa. 2008. An algorithm for non-deterministic object distribution in P systems and its implementation in hardware. In Proceedings of the 9th International Workshop on Membrane Computing (WMC’08) (LNCS), D. W. Corne, P. Frisco, Gh. Păun, G. Rozenberg, and A. Salomaa (Eds.), Vol. 5391. Springer, 325–354. [72] Van Nguyen, David Kearney, and Gianpaolo Gioiosa. 2008. An implementation of membrane computing using reconfigurable hardware. Comput. Inform. 27, 3 (2008), 551–569. [73] Van Nguyen, David Kearney, and Gianpaolo Gioiosa. 2009. A region-oriented hardware implementation for membrane computing applications. In Proceedings of the 10th International Workshop on Membrane Computing (WMC’09) (LNCS), Gh. Păun, M. J. Pérez-Jiménez, A. Riscos-Núñez, G. Rozenberg, and A. Salomaa (Eds.), Vol. 5957. Springer, 385–409. [74] Van Nguyen, David Kearney, and Gianpaolo Gioiosa. 2010. An extensible, maintainable and elegant approach to hardware source code generation in Reconfig-P. J. Logic Algeb. Program. 79, 6 (2010), 383–396. [75] Taishin Y. Nishida. 2006. A Membrane Computing Model of Photosynthesis. Springer, 181–202. [76] David Orellana-Martín, Miguel A. Martínez-del-Amor, Luis Valencia-Cabrera, Bosheng Song, Linqiang Pan, and Mario J. Pérez-Jiménez. 2020. P systems with symport/antiport rules: When do the surroundings matter? Theoret. Comput. Sci. 805 (2020), 206–217. [77] David Orellana-Martín, Luis Valencia-Cabrera, Agustín Riscos-Núñez, and Mario J. Pérez-Jiménez. 2019. Minimal cooperation as a way to achieve the efficiency in cell-like membrane systems. J. Memb. Comput. 1, 2 (2019), 85–92. [78] Ana Pavel, Octavian Arsene, and Catalin Buiu. 2010. Enzymatic numerical P systems—A new class of membrane computing systems. In Proceedings of the 5th International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA’10), A. K. Nagar, R. Thamburaj, K. Li, Z. Tang, and R. Li (Eds.). IEEE, 1331–1336. [79] Ignacio Pérez-Hurtado, Miguel Á. Martínez-del-Amor, Gexiang Zhang, Ferrante Neri, and Mario J. Pérez-Jiménez. 2020. A membrane parallel rapidly-exploring random tree algorithm for robotic motion planning. Integ. Comput.- Aided Eng. 27 (2020), 121–138. [80] Ignacio Pérez-Hurtado, David Orellana-Martín, Gexiang Zhang, and Mario J. Pérez-Jiménez. 2019. P-Lingua in two steps: Flexibility and efficiency. J. Memb. Comput. 1, 2 (01 June 2019), 93–102. [81] Ignacio Pérez-Hurtado, Mario J. Pérez-Jiménez, Gexiang Zhang, and David Orellana-Martín. 2018. Simulation of rapidly-exploring random trees in membrane computing with P-lingua and automatic programming. Int. J. Comput. Commun. Contr. 13, 6 (2018), 1007–1031. [82] Biljana Petreska and Christof Teuscher. 2003. A reconfigurable hardware membrane system. In Proceedings of the International Workshop on Membrane Computing (WMC’03) (Lecture Notes in Computer Science), Carlos Martín-Vide, Giancarlo Mauri, Gheorghe Paun, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 2933. Springer, 269–285. [83] Andrew Pohorille and David Deamer. 2009. Self-assembly and function of primitive cell membranes. Res. Microbiol. 160, 7 (2009), 449–456. [84] Gheorghe Păun. 2000. Computing with membranes. J. Comput. Syst. Sci. 61, 1 (2000), 108–143. [85] Gheorghe Păun. 2002. Membrane Computing: An Introduction. Springer-Verlag, Berlin. [86] Gheorghe Păun and Radu A. Păun. 2006. Membrane computing and economics: Numerical P systems. Fundam. Inform. 73, 1–2 (2006), 213–227. [87] Gheorghe Păun, Grzegorz Rozenberg, and Arto Salomaa (Eds.). 2009. The Oxford Handbook of Membrane Computing. Oxford University Press. [88] Juan Quiros, Sergey Verlan, Julian Viejo, Alejandro Millán, and Manuel J. Bellido. 2016. Fast hardware implementations of static P systems. Comput. Inform. 35, 3 (2016), 687–718. [89] Raúl Reina-Molina, Daniel Díaz-Pernil, and Miguel A. Gutiérrez-Naranjo. 2011. Integer linear programming for tissue-like P systems. In Proceedings of the 9th Brainstorming Week on Membrane Computing,M.A.Martínez-delAmor, Gh. Păun, I. Pérez-Hurtado, F. J. Romero-Campero, and L. Valencia-Cabrera (Eds.). [90] Haina Rong, Kang Yi, Gexiang Zhang, Jianping Dong, Prithwineel Paul, and Zhiwei Huang. 2019. Automatic implementation of fuzzy reasoning spiking neural P systems for diagnosing faults in complex power systems. Complexity 2019 (2019), 2635714:1–2635714:16. [91] Grzegorz Rozenberg and Arto Salomaa (Eds.). 1997. Handbook of Formal Languages.Vol.1–3.Springer. [92] Eduardo Sánchez-Karhunen and Luis Valencia-Cabrera. 2019. Modelling complex market interactions using PDP systems. J. Memb. Comput. 1, 1 (2019), 40–51. [93] James E. Smith. 1984. Decoupled access/execute computer architectures. ACM Trans. Comput. Syst. 2, 4 (1984), 289– 308. [94] Yasuhiro Suzuki, Yoshi Fujiwara, Junji Takabayashi, and Hiroshi Tanaka. 2000. Artificial life applications of a class of P systems: Abstract rewriting systems on multisets. In Proceedings of the Workshop on Membrane Computing - Multiset Processing (LNCS), C. S. Calude, Gh. Păun, G. Rozenberg, and A. Salomaa (Eds.). Springer, 299–346. [95] Yasuhiro Suzuki and Hiroshi Tanaka. 2006. Modeling p53 Signaling Pathways by Using Multiset Processing. Springer Berlin, 203–214. [96] El-Ghazali Talbi, Sanaz Mostaghim, Tatsuya Okabe, Hisao Ishibuchi, Günter Rudolph, and Carlos A. Coello Coello. 2008. Parallel approaches for multiobjective optimization. In Multiobjective Optimization, Interactive and Evolutionary Approaches (LNCS), J. Branke, K. Deb, K. Miettinen, and R. Slowinski (Eds.), Vol. 5252. Springer, 349– 372. [97] El-Ghazali Talbi. 2018. A unified view of parallel multi-objective evolutionary algorithms. J. Parallel Distrib. Comput. 133 (2018), 349–358. [98] Jorge A. Tejedor, Luis Fernández, Fernando Arroyo, and Sandra Gómez-Canaval. 2007. Algorithm of rules applications based on competitiveness of evolution rules. In Proceedings of the 8th Workshop on Membrane Computing,G. Ekeftherakis, P. Kefalas, and Gh. Păun (Eds.). 567–580. [99] Jorge A. Tejedor, Luis Fernández, Fernando Arroyo, and Abraham Gutiérrez. 2007. Algorithm of active rule elimination for application of evolution rules. In Proceedings of the 8th WSEAS International Conference on Evolutionary Computing, Akshai Aggarwal (Ed.). 259–267. [100] Jorge A. Tejedor, Abraham Gutiérrez, Luis Fernández, Fernando Arroyo, Ginés Bravo, and Sandra Gómez-Canaval. 2007. Optimizing evolution rules application and communication times in membrane systems implementation. In Proceedings of the 8th International Workshop on Membrane Computing (WMC’07) (Lecture Notes in Computer Science), George Eleftherakis, Petros Kefalas, Gheorghe Păun, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 4860. Springer, 298–319. [101] Roman Trobec, Marián Vajteršic, and Peter Zinterhof. 2009. Parallel Computing. Springer. [102] Christos Tsotskas, Timoleon Kipouros, and Anthony Mark Savill. 2014. The design and implementation of a GPUenabled multi-objective Tabu-Search intended for real world and high-dimensional applications. Proced. Comput. Sci. 29 (2014), 2152–2161. [103] Ganesh Venkatesh, Jack Sampson, Nathan Goulding, Sravanthi Kota Venkata, Michael Bedford Taylor, and Steven Swanson. 2011. QSCORES: Trading dark silicon for scalable energy efficiency with quasi-specific cores. In Proceedings of the 44th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO’11). 163–174. [104] Sergey Verlan. 2010. Study of Language-Theoretic Computational Paradigms Inspired by Biology. Habilitation thesis, Université Paris Est. [105] Sergey Verlan. 2013. Using the formal framework for P systems. In Proceedings of the 14th International Conference on Membrane Computing (CMC’13) (Lecture Notes in Computer Science), Artiom Alhazov, Svetlana Cojocaru, Marian Gheorghe, Yurii Rogozhin, Grzegorz Rozenberg, and Arto Salomaa (Eds.), Vol. 8340. Springer, 56–79. [106] Sergey Verlan and Juan Quiros. 2012. Fast hardware implementations of P systems. In Proceedings of the 13th International Conference on Membrane Computing (CMC’12) (Lecture Notes in Computer Science), Erzsébet Csuhaj-Varjú, Marian Gheorghe, Grzegorz Rozenberg, Arto Salomaa, and György Vaszil (Eds.), Vol. 7762. Springer, 404–423. [107] Tao Wang, Gexiang Zhang, and Mario J. Pérez-Jiménez. 2015. Fuzzy membrane computing: Theory and applications. Int. J. Comput. Commun. Contr. 10 (2015), 904–935. [108] Tao Wang, Gexiang Zhang, Junbo Zhao, Zhenyou He, Jun Wang, and Mario J. Pérez-Jiménez. 2015. Fault diagnosis of electric power systems based on fuzzy reasoning spiking neural P systems. IEEE Trans. Power Syst. 30, 3 (2015), 1182–1194. [109] Xueyuan Wang, Gexiang Zhang, Ferrante Neri, Tao Jiang, Junbo Zhao, Marian Gheorghe, Florentin Ipate, and Raluca Lefticaru. 2016. Design and implementation of membrane controllers for trajectory tracking of nonholonomic wheeled mobile robots. Integ. Comput.-Aided Eng. 23, 1 (2016), 15–30. [110] Xueyuan Wang, Gexiang Zhang, Junbo Zhao, Haina Rong, Florentin Ipate, and Raluca Lefticaru. 2015. A modified membrane-inspired algorithm based on particle swarm optimization for mobile robot path planning. Int. J. Comput. Commun. Contr. 10, 5 (2015), 732–745. [111] Nicholas Wilt (Ed.). 2013. The CUDA Handbook: A Comprehensive Guide to GPU Programming. Addison Wesley. [112] Erica Wiseman. 2016. Next generation computing. National Research Council of Canada/Gov. of Canada (2016). https: //cradpdf.drdc-rddc.gc.ca/PDFS/unc268/p805200_A1b.pdf. [113] Zihan Xu, Matteo Cavaliere, Pei An, Sarma Vrudhula, and Yu Cao. 2014. The stochastic loss of spikes in spiking neural P systems: Design and implementation of reliable arithmetic circuits. Fund. Inform. 134, 1–2 (2014), 183–200. [114] Jianying Yuan, Dequan Guo, Gexiang Zhang, Prithwineel Paul, Ming Zhu, and Qiang Yang. 2019. A resolution-free parallel algorithm for image edge detection within the framework of enzymatic numerical P systems. Molecules 24, 7 (2019). [115] Xiangxiang Zeng, Henry Adorna, Miguel A. Martínez-del-Amor, Linqiang Pan, and Mario J. Pérez-Jiménez. 2010. Matrix representation of spiking neural P systems. In Proceedings of the 11th International Conference on Membrane Computing (CMC’10) (LNCS), Marian Gheorghe, Thomas Hinze, Gheorghe Păun, Grzegorz Rozenberg, and Arto Salomaa (Eds.). 377–391. [116] Gexiang Zhang, Jixiang Cheng, Marian Gheorghe, and Qi Meng. 2013. A hybrid approach based on differential evolution and tissue membrane systems for solving constrained manufacturing parameter optimization problems. Appl. Soft Comput. 13, 3 (2013), 1528–1542. [117] Gexiang Zhang, Marian Gheorghe, Linqiang Pan, and Mario J. Pérez-Jiménez. 2014. Evolutionary membrane computing: A comprehensive survey and new results. Inform. Sci. 279 (2014), 528–551. [118] Gexiang Zhang, Mario J. Pérez-Jiménez, and Marian Gheorghe. 2017. Real-life Applications with Membrane Computing (1st ed.). Springer Publishing Company, Incorporated. [119] Gexiang Zhang, Haina Rong, Ferrante Neri, and Mario J. Pérez-Jiménez. 2014. An optimization spiking neural P system for approximately solving combinatorial optimization problems. Int. J. Neural Syst. 24, 5 (2014), 1440006. [120] Weihang Zhu, Ashraf Yaseen, and Yaohang Li. 2011. DEMCMC-GPU: An efficient multi-objective optimization method with GPU acceleration on the Fermi architecture. New Gen. Comput. 29, 2 (2011), 163–184.