Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations

139 Fundamenta Informaticae 164 (2019) 139–155 DOI 10.3233/FI-2019-1759 Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations Jym Paul Carandang, Francis George C. Cabarle*, Henry Natividad Adorna, Nestine Hope S. Hernandez Algorithms and Complexity Department of Computer Science University of the Philippines Diliman, Philippines [email protected], {fccabarle, hnadorna, nshernandez}@up.edu.ph Miguel ´ Angel Mart´ ınez-del-Amor Research Group on Natural Computing Department of Computer Science and Artificial Intelligence Universidad de Sevilla, Spain [email protected] Abstract. Spiking Neural P system is a computing model inspired on how the neurons in a living being are interconnected and exchange information. As a model in embrane computing, it is a non-deterministic and massively-parallel system. The latter makes GPU a good candidate for accelerating the simulation of these models. A matrix representation for systems with and without delay have been previously designed, and algorithms for simulating them with deterministic systems was also developed. So far, non-determinism has been problematic for the design of parallel simulators. In this work, an algorithm for simulating non-deterministic spiking neural P system with delays is presented. In order to study how the simulations get accelerated on a GPU, this algorithm was implemented in CUDA and used to simulate non-uniform and uniform solutions to the Subset Sum problem as a case study. The analysis is completed with a comparison of time and space resources in the GPU of such simulations. *Address for correspondence: Algorithms and Complexity, Department of Computer Science, University of the Philippines Diliman, Diliman 1101 Quezon City, Philippines 140 J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations Keywords: Membrane Computing; Spiking Neural P systems; Matrix Representation; CUDA; GPU; Subset Sum 1. Introduction Spiking neural P system (or SN P system) is a computing model under membrane computing [1] that is inspired on neurons in a living being where each neuron sends electric signals to each other to deliver information [2]. The model consists of computing units called neurons which contain spikes. The neurons are interconnected via synapses and pass spikes to each other based on a set of rules. The model works in a parallel and non-deterministic manner, i.e., each neuron “processes” their spikes simultaneously and non-deterministically chooses which rule to apply in case multiple rules can be selected. Simulators for SN P systems, as for other types of P systems, will allow to validate the models experimentally, as well as to analyze how well they perform in current technologies within High Performance Computing. SN P systems have been investigated in many aspects, such as their computing power (under various modes or context) in [2, 3, 4, 5], their computing efficiency in [6, 7, 8] to name a few. Further, many variants of SN P systems have been introduced based on many functions and structures of the biological neuron, such as [9, 10, 11, 12, 13] and more recently [14, 15, 16] to name a few. There are already efforts in creating simulators for SN P systems like [17, 18, 19, 20] with earlier versions of the CuSNP simulator in [21], but they are only designed for deterministic systems. Non-determinism is a feature that usually has a strong impact in the simulators. So far there have been efforts on reducing this cost, as the DND algorithm for hardware implementation of transition P systems [22]. In [23], it is discussed that the DCBA algorithm, designed for PDP systems, includes a phase that handles non-determinism. This phase is inherently sequential and breaks parallelism when offloaded to the GPU. In this paper, an algorithm for managing non-determinism in the simulation of SN P systems is presented. This is mainly done by extending representations in the previous work from vectors to matrices where each row corresponds to a different configuration of the system. The algorithm was then implemented as a case study and used to simulate uniform and non-uniform solutions to the Subset Sum problem as presented in [7]. Time and space usage for both simulations are then compared and analyzed to give insights regarding how to further accelerate simulations of non-deterministic spiking neural P systems. An algorithm for generating Spiking Vectors is also presented. The algorithm works by first finding all possible rules that can fire at that time step. However, not all rules can fire so different valid spiking vectors that follows the semantics of an SN P system are generated and collected as rows of a Spiking Matrix. The paper is structured as follows: Section 2 provides a formal definition of the spiking neural P system and some matrix representation is presented; Section 2.2 presents CUDA programming paradigm; Section 3 presents the algorithm for simulating non deterministic SN P systems as well as generating spiking vectors as well as matrix representation used in this report; Section 4 shows the results of implementing the algorithm and simulating different systems; and Section 5 provides conclusion as well as some future work. J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations 141 2. Preliminaries 2.1. Spiking neural P systems We assume that the reader has basic knowledge in formal language and automata theory, and membrane computing. We only briefly mention notations and definitions relevant to what follows. More information can be found in various monographs, e.g. [24] and [25]. Let Vbe an alphabet, V*is the set of all finite strings over V, the empty string is denoted as λ, and V+=V∗− {λ}.If ais a symbol in V, then a0=λ. A language defined by a regular expression Eover Vis denoted as L(E). Only the following are regular expressions over V: (i) the empty string λ, (ii) each a∈V, and (iii)E1E2– concatentation – , E1∪E2– union –, E+ 1– positive closure– where E1and E2are regular expressions over V. We extend this definition to languages by regular expressions E1and E2 in a natural way, i.e., L(λ) = {λ},L(a) = {a}for a∈V, and L(E1E2) = L(E1)L(E2), L(E1∪E2) = L(E1)∪L(E2), L(E+ 1) = E+ 1. If V={a}, i.e., singleton set, we write a∗and a+instead of {a}∗and {a+}, respectively, for simplicity. Next, the model of spiking neural P systems from [26] is briefly covered. Spiking neural P system Πof degree m≥1is of the form: Π=(O, σ1, ..., σm, syn, in, out) 1. O={a}is the alphabet containing a single symbol (the spike); 2. σ1, ..., σmare neurons, of the form σj= (αj, Rj),1≤j≤m, where: (a) αj≥0is the initial number of spikes contained in σj. (b) Rjis a finite set of rules of the following two forms: (i) (Spiking Rule) E/ac→ap;dwhere Eis a regular expression over Oand c≥p≥ 1, d ≥0. (ii) (Forgetting Rule) as→λ, for s≥1,with the restriction that for each rule E/ac→ ap;dof type (i) from Rj,as6∈ L(E); 3. syn ⊆ {1,2, ..., m}×{1,2, ..., m}with i6=jfor all (i, j)∈syn, 1≤i, j ≤m(synapses between neurons); 4. in, out ∈ {1,2, ..., m}indicate the input and the output neurons, respectively. An SN P system whose spiking rules have p= 1 is said to be of the standard type (non-extended). A spiking rule is applied as follows: if a neuron σjcontains kspikes, and ak∈L(E), k ≥c, then the rule r≡E/ac→ap;d, for r∈Rj, can be applied. This means cspikes are removed so that k−c spikes remain in σj, the neuron is then fired and produces pspikes after dtime steps. That is, spikes are fired on step t+d, where tis the current time step of the computation. If d= 0, the spikes are fired immediately. From step tto t+dneuron σhas not fired the spike yet and is said to be closed, i.e., σjcannot receive spikes from other neurons connected to it. If neurons with a synapse to σjfire, the spikes are lost. At step t+d,σjis open again to receive spikes, and at t+d+ 1,σjcan apply a rule. A forgetting rule is applied as follows: if σjcontains exactly sspikes, then the rule r≡as→λ from Rjcan be applied, meaning all of the sspikes are removed from σj. 142 J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations Rules of type (i) where E=accan be written in the shortened form of ac→ap;d. In the case two or more rules of a neuron are applicable at the same step, only one rule is applied and is non-deterministically chosen. A configuration of the system at step tis denoted as Ct= hc1/k1, . . . , cm/kmi, where σjcontains cj≥0spikes and remains closed for kjmore steps, for 1≤j≤m. The initial configuration of the system is therefore C0=hc1/0, . . . , cm/0i. Rule application provides us a transition from one configuration to another. A computation is any (finite or infinite) sequence of configurations such that: (a) the first term is the initial configuration C0; (b) for each n≥1, the nth configuration of the sequence is obtained from the previous configuration in one transition step; and (c) if the sequence is finite (called halting computation) then the last term is a halting configuration, i.e., a configuration where all neurons are open and no rule can be applied. Two common ways to obtain the output of an SN P system are as follows: (1) obtaining the time interval between exactly the first two steps when the output neuron σout spikes, e.g. number n=t2−t1 is computed, where σout produced its first two spikes at steps t1and t2; (2) counting the number of spikes produced by σout until the system halts. In this work, systems that produce their output as given in (2) are considered, as done in [7]. In this paper, the focus is on SN P systems solving the Subset Sum problem mainly, and how the uniform and non-uniform solution to this problem compares to each other. When a system is constructed in such a way that it depends directly at an instance of a problem, it is said that the solution is non-uniform, otherwise it is uniform. For the Subset Sum problem for example, it suffices to know the number of integers in the set for the uniform solution while it is needed to know both the number of integers in the set as well as each integer for the non-uniform solution as the number of neurons will also depend on this [7]. 2.2. Parallel computing with CUDA GPUs Today Graphics Processing Units (GPUs) feature thousands of processing cores that can be harnessed for parallel computing. In this concern, the Compute Unified Device Architecture (in short, CUDA) is a parallel programming computing platform and application programming interface model developed by NVIDIA [27] [28]. This technology allows software developers to use CUDA-enabled GPUs for scientific computing, also known as general-purpose computations on the GPU (in short, GPGPU), or simply, GPU computing. Functions that execute in the GPU, known as kernel functions, are executed by one or more threads arranged in thread blocks. The thread blocks are then arranged in a grid of thread blocks. It is the task of the developer to arrange threads in such a thread hierarchy, i.e., in a block of threads, and a grid of thread blocks. In the CUDA programming model, the GPU is often referred to as the device, while the CPU is referred to as the host. The host performs the kernel function calls to be executed on the device. CUDA uses a single thread, multiple data (in short, STMD) paradigm. In this paradigm, threads execute similar code, while allowing such threads to access multiple (possibly different values of) data. Aside from the thread hierarchy, CUDA also implements a memory hierarchy similar to how there exist memory hierarchies in the CPU. For more details of CUDA thread and memory hierarchy can be found in [27] [28]. J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations 143 Linear algebra operations fit very well with GPU computing paradigm. In fact, there are several optimized libraries for that freely available. Indeed, GPUs can handle very efficiently matrix and vector operations, since their design is adapted for that because of the nature of graphics operations. 2.3. Related work on simulating P systems with GPUs Parallel simulation of SN P systems has been based so far on the matrix representation introduced originally in [29]. The aim of it is to express the syntactic elements in form of matrices, so that the semantics are mainly represented by matrix operations. Given that GPUs are good at this kind of operations, parallel implementations are expected to accelerate the sequential counterpart. This representation was later extended to admit rules with delays in [18]. The framework of parallel simulators for SN P systems is called CuSNP [17], and has been also iteratively improved to both admit more ingredients as delays [18] and more types of regular expressions. These simulators work on Python and C++, with the interplay of CUDA. GPU simulators for non-neural-like P systems are also available. Specifically, first GPU simulators were available for cell–like P systems (with active membranes) [30], and these were followed by other parallel simulators for tissue–like P systems (with cell division), enzymatic numerical P systems and Population Dynamics P (PDP) systems [31]. Non-determinism was taken into account when designing them. For example, the simulator for P systems with active membranes assumes that the simulated models are confluent; since all computations halt and give the same answer, the simulator can choose only one path on the fly using lightweight strategies (for example, giving more priorities to evolution rules than to division rules) [30]. The strategy to handle non-determinism in PDP systems was not to get rid of it in the simulator, but to implement strategies directly in the simulation algorithm DCBA [23]. Although this model includes probabilities in the application of rules, this is done by computing multinomial random variates for rules having the very same left-hand side. Thus, non-determinism takes place for selecting rules having overlapping left-hand sides. For this purpose, the algorithm includes a specific phase to take control on this, by imposing a random order of remaining rules to be applied, and looping over them by assigning the remaining objects. This phase was found to be a bottleneck on the stress tests carried out on GPUs [23]. 3. Algorithms for CuSNP In this section, the modifications to the matrix representation given in [29], later extended for delays in [18], are first introduced. They aim at allowing the simulation of SN P systems with non-determinism. After the modifications are presented, the algorithms for the simulation in CuSNP are discussed. 3.1. Matrix representation In order to show in a illustrative way the following matrix representation, the SN P system Πepresented in Figure 1 is used as an example. Πecan be formally defined as follows: Πe= ({a}, σ1, σ2, syn) 144 J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations a2 (1) a2 → a2 (2) a2 → a a (3) a → a (4) a → a;1 12 Figure 1: Πe, a spiking neural P system with 2 neurons where 1. σ1= (2,{a2/a2→a2; 0, a2/a2→a; 0}), 2. σ2= (1,{a/a →a; 0, a/a →a; 1}), 3. syn ={(1,2),(2,1)} The following definitions are used to construct the matrix representation. Let Πbe an SN P system with delays having mneurons and nrules in total. Note that nis calculated as n=Pm j=1 |Rj|, where σj= (αj, Rj)are the neurons of Π. Let assume a total order of the rules in Π, so that each rule can be unequivocally identified with a number 1≤i≤n, as ri. Definition 3.1. The Configuration Vector C(t)=hct 1, ..., ct mi, where ct jis the amount of spikes in σjat time t, as done in [29]. Definition 3.2. The tth Status Vector is denoted by St(t)=hstt 1, ..., stt mifrom [18] such that, stt j=(1,if neuron σjis open in step t, 0,if neuron σjis closed in step t. Definition 3.3. The tth Countdown Vector is represented by D(t)=hdt 1, ..., dt nifrom [18] such that, dt i=     −1,if rule riis not fired at step t, 0,if spikes of rule riare released at step t, ≥1,if rule riis fired but has a delay at step t Note that ditells us the time remaining before rule rireleases its spikes. Definition 3.4. The System Configuration Conf(t)at time step tis finally defined by the following: Conf(t)= (C(t), St(t), D(t)) For the example Πe, with degree m= 2, the initial Configuration Vector is C(0) =h2,1i. The initial Status Vector is St(0) =h1,1isince all the neurons are initially open. The initial countdown vector D(0) =h−1,−1,−1,−1isince all of the rules are not fired yet. Thus, Conf(0) = (C(0), St(0), D(0)) is the initial System Configuration. J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations 145 Definition 3.5. The Spiking Matrix of a system Πat time step tcan be defined as follows: S(t)=hs(t) ki iqt×n Where: s(t) ki =(1,if for rule riin Rjof σj, E is satisfied, dt i=−1, stt j= 1, and it is applied in instance k 0,otherwise. That is to say that each row of the Spiking Matrix Sis composed of individual valid spiking vectors as defined in [18]. In our case, qtis the number of valid spiking vectors from a single configuration or instance Conf(t). Definition 3.6. An Indicator Matrix of a system Πat time step tcan be defined as follows: IM(t)=him(t) ki iqt×n where im(t) ki =(1,if rule riwith delay dfired at time step s, where s+d=t 0,otherwise. Like the Spiking Matrix S, the Indicator Vector is also extended into a matrix where each row represents a single indicator vector. In other words, the value is 1 if the corresponding rule fired and is done with its delay. Definition 3.7. AStatus Matrix of a system Πat time step tcan be defined as follows: SM(t)=hsm(t) kj iqt×m where sm(t) kj =(1,if neuron σjis open at instance kin time step t 0,otherwise. Definition 3.8. The Transition Matrix from [29] of a system Πcan be defined as follows: TM =htmijin×m where tmij =     p, if rule riis in neuron σs(s6=j and (s, j)∈syn)and produces pspikes −c, if rule riis in neuron σjand consumes cspikes 0,otherwise. 146 J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations Definition 3.9. From Definition 3.8 for TM, two more matrices, TM+and TM−, are defined as follows TM+=htm+ ijin×mand TM−=htm− ijin×m where tm+ ij =(tmij,if tmij >0 0,otherwise. tm− ij =(tmij,if tmij <0 0,otherwise. That is to say that TM+contains the positive values in TM (i.e., the number of spikes produced by the rules) and TM−contains the negative values (i.e., the number of spikes consumed) Definition 3.10. The Net Gain Matrix at a time step tcan be defined as: NG(t)=hC(t+1) w−C(t)iqt Each row w, 1≤w≤qt, is a single Net Gain Vector [18]. The idea for NG(t)is to store the net effect of spikes consumed and produced for each neuron at step k. Unlike in [29] where delays in rules are not represented, here we need to check for spikes produced by neurons that become open at step k, hence multiplying IM(t)to TM+ n×m. For the above definition, NG(t)can be computed by: NG(t) qt×m=SM(t) qt×m⊗(IM(t) qt×n∗TM+ n×m) + S(t) qt×n∗TM− n×m where ⊗is the elementwise multiplication for matrix. Continuing with the example, the following are the matrices for Πein Figure 1: TM+=      0 2 0 1 1 0 1 0       TM−=      −2 0 −2 0 0−1 0−1       S(1) =     1010 1001 0110 0101      IM(1) =     1010 1000 0110 0100      SM(1) =     1 1 1 0 1 1 1 0      NG(1) =     1 1 1 0 1 1 1 0      ⊗          1010 1000 0110 0100      ∗     0 2 0 1 1 0 1 0           +          1010 1001 0110 0101      ∗     −2 0 −2 0 0−1 0−1           =     −1 1 −2−1 −1 0 −2−1      C(1) =C(0) +NG(1) J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations 147 C(1) =     2 1 2 1 2 1 2 1      +     −1 1 −2−1 −1 0 −2−1      =     1 2 0 0 1 1 0 0      3.2. Main simulation algorithm Next, the algorithms that make use of the definitions provided above are provided. Algorithm 1 refers to the here-presented algorithm for simulating non-deterministic spiking neural P systems, and Algorithm 2 provides the main loop. Algorithm 1 Core simulation algorithm for non-deterministic SN P systems Require: Queue Q, containing “un-explored” configurations 1: Get a system configuration Conf(t)from Q 2: Get the Countdown vector D(t)and Status Matrix St(t)from Conf(t) 3: Generate the Spiking Matrix S(t) qt×n, where qtis the number of possible spiking vectors. 4: Generate the Indicator Matrix IM(t) qt×n 5: NG(t) q×m←St(t)⊗(IM(t)∗TM+) + S(t)∗T M− 6: D(t+1) ←D(t) 7: for each dt+1 i∈D(t+1) do 8: dt+1 i←MAX(dt+1 i−1,−1) 9: end for 10: for k←1. . . qtdo 11: St(t+1) k←Update the Status Vector according to [S(t) k]nand St(t) 12: D(t+1) k←Update the Countdown Vector according to [S(t) k]nand D(t+1) 13: Generate a new system configuration Conf(t+1) k= (C(t)+NG(t) k, St(t+1) k, D(t+1) k). 14: if Conf(t+1) kis not a halting configuration then 15: Add Conf(t+1) kto Q 16: end if 17: end for 18: return Q The key of this algorithm is the use of a queue (i.e., a first-in, first out data structure implemented with a 1D array) to store the configurations generated in Algorithm 1, so the computation tree is processed in a Breadth-First manner. These steps are repeated in order to simulate the system to a definite number of iterations or until the computation halts and no configuration remains in the queue. The Indicator Matrix, Countdown Matrix, and Status Matrix are generated based on their definitions above using the Spiking Vector. The next section will present how the spiking matrix is created. In particular, non-determinism of SN P systems is simulated in Algorithm 1 when generating the spiking 154 J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations [12] Song T, Pan L, P˘ aun G. Spiking neural P systems with rules on synapses. Theoretical Computer Science. 2014;529:82–95. URL https://doi.org/10.1016/j.tcs.2014.01.001. [13] Cabarle FGC, Adorna HN, P´ erez-Jim´ enez MJ, Song T. Spiking neural P systems with structural plasticity. Neural Computing and Applications. 2015;26(8):1905–1917. doi:10.1007/s00521-015-1857-4. [14] Song T, Pan L. Spiking neural P systems with request rules. Neurocomputing. 2016;193:193–200. URL https://doi.org/10.1016/j.neucom.2016.02.023. [15] Cabarle FGC, Adorna HN, Jiang M, Zeng X. Spiking neural P systems with scheduled synapses. IEEE Trans Nanobiosci. 2017;16(8):792–801. doi:10.1109/TNB.2017.2762580. [16] Wu T, P˘ aun A, Zhang Z, Pan L. Spiking neural P systems with polarizations. IEEE Trans Neural Net Learning Sys. 2017;. [17] Cabarle FGC, Adorna HN, Mart´ ınez-del-Amor M ´ A, P´ erez-Jim´ enez MJ. Improving GPU Simulations of Spiking Neural P Systems. Romanian Journal of Information Science and Technology. 2012;15(1):5–20. URL http://www.imt.ro/romjist/Volum15/Number15_1/cuprins15_1.htm. [18] Carandang JPA, Villaflores JMB, Cabarle FGC, Adorna HN, Mart´ ınez-del-Amor M ´ A. CuSNP: Spiking Neural P Systems Simulators in CUDA. Romanian Journal of Information Science and Technology. 2017;20(1):57–70. URL http://2016.asiancmc.org/. [19] Carandang JPA, Villaflores JMB, Cabarle FGC, Adorna H, Mart´ ınez-del-Amor M ´ A. CuSNP: Spiking Neural P Systems Simulators in CUDA. 5th Asian Conference on Membrane Computing. 2016. [20] Carandang JPA, Villaflores JMB, Cabarle FGC, Adorna H, Mart´ ınez-del-Amor M ´ A. Improving simulations of Spiking Neural P Systems in NVIDIA CUDA GPUs: CuSNP. In: 14th Brainstorming Week on Membrane Computing. Seville, Spain: F´ enix Editora; 2016. pp. 135–150. URL http://www.gcn.us. es/14bwmc_proceedings. [21] Carandang JPA, Villaflores JMB, Cabarle FGC, Adorna H, Mart´ ınez-del-Amor M ´ A. CuSNP version 06.06.16; 2016. http://aclab.dcs.upd.edu.ph/productions/software/cusnp_v060616. [22] Nguyen V, Kearney D, Gioiosa G. An Algorithm for Non-deterministic Object Distribution in P Systems and Its Implementation in Hardware. In: Corne DW, Frisco P, P˘ aun G, Rozenberg G, Salomaa A, editors. Membrane Computing. Berlin, Heidelberg: Springer Berlin Heidelberg; 2009. pp. 325–354. doi:10.1007/978-3-540-95885-7 24. [23] Mart´ ınez-Del-Amor, Miguel ´ Angel and Mac´ ıas-Ramos, Luis Felipe and Valencia-Cabrera, Luis and P´ erez- Jim´ enez, Mario J. Parallel Simulation of PDP Systems: Updates And Roadmaps. Natural Computing 2016; 15(4):565–573. doi:10.1007/s11047-016-9566-1 [24] P˘ aun G. Membrane Computing: An Introduction. Springer-Verlag Berlin Heidelberg: Springer Berlin Heidelberg; 2002. doi:10.1007/978-3-642-56196-2. [25] P˘ aun G, Rozenberg G, Salomaa A, editors. The Oxford Handbook of Membrane Computing. Oxford Univeristy Press; 2010. ISBN:0199556679, 9780199556670. [26] Ionescu M, P˘ aun G, Yokomori T. Spiking Neural P Systems. Fundam Inf. 2006;71(2-3):279–308. [27] NVIDIA. CUDA C Programming Guide; 2015. Accessed: 2015-11-19. http://docs.nvidia.com/ cuda/cuda-c-programming-guide. [28] Kirk DB, Hwu WmW. Programming Massively Parallel Processors: A Hands-on Approach. 2nd ed. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.; 2013. ISBN:9780123914187, 9780124159921. J.P. Carandang et al. /Handling Non-determinism in Spiking Neural P Systems: Algorithms and Simulations 155 [29] Zeng X, Adorna H, Mart´ ınez-del-Amor M ´ A, Pan L, P´ erez-Jim´ enez MJ. In: Gheorghe M, Hinze T, P˘ aun G, Rozenberg G, Salomaa A, editors. Matrix Representation of Spiking Neural P Systems. Berlin, Heidelberg: Springer Berlin Heidelberg; 2011. pp. 377–391. doi:10.1007/978-3-642-18123-8 29. [30] Cecilia JM, Garc´ ıa JM, Guerrero GD, Mart´ ınez-del-Amor M ´ A, P´ erez-Hurtado I, P´ erez-Jim´ enez MJ. Simulation of P systems with active membranes on CUDA. Briefings in Bioinformatics. 2010;11(3):313–322. doi:10.1093/bib/bbp064. [31] Mart´ ınez-del-Amor M ´ A, Garc´ ıa-Quismondo M, Mac´ ıas-Ramos LF, Valencia-Cabrera L, Riscos-N´ u˜ nez A, P´ erez-Jim´ enez MJ. Simulating P Systems on GPU Devices: A Survey. Fundamenta Informaticae. 2015;136:269–284. doi:10.3233/FI-2015-1157. [32] Mart´ ınez-del-Amor M ´ A, Orellana-Mart´ ın D, Cabarle FGC, P´ erez-Jim´ enez MJ, Adorna HN. Sparse-matrix Representation of Spiking Neural P Systems for GPU. In: Proc. of 15th Brainstorming Week on Membrane Computing. Seville, Spain: F´ enix Editora; 2017. pp. 161–170. URL http://www.gcn.us.es/15bwmc_ proceedings.