BACKGROUND
systems biology is a multidisciplinary research field relying on the cross-talk between mathematical, computational and experimental tools to investigate the functioning of complex biological systems, and to predict how they might behave in both physiological and perturbed conditions.
to this aim, different computational methods—e.g., parameter estimation, sensitivity analysis or reverse engineering  <cit> —are usually exploited to define or calibrate the mathematical model that describes the system of interest.
these methods require the execution of a large number of simulations, each one generally corresponding to a distinct model structure or parameterization, that is, to a different set of molecular interactions or to different initializations of the species amounts and/or reaction constants.
as a result, the computational burden required by these computational analyses can rapidly overtake the capabilities of central processing units , therefore limiting in-depth computational investigations to small-scale models consisting in a few tens of reactions and molecular species at most.
general-purpose graphics processing units  can be exploited to overcome these drawbacks.
indeed, they are parallel multi-core co-processors that are drawing an ever-growing attention by the scientific community, since they give access to tera-scale performances on common workstations .
as such, they can markedly decrease the running times required by traditional cpu-based software, still maintaining low-costs and energetic efficiency.
as a matter of fact, in the latter years gpus have been widely adopted as an alternative approach to classic parallel architectures for the parallelization of computational methods in systems biology, computational biology and bioinformatics  <cit> .
in this work we propose lassie , a novel gpu-accelerated software designed to simulate large-scale reaction-based models of cellular processes, consisting in hundreds or thousands of reactions and molecular species.
an example of killer-application of lassie would consist in the simulation of rule-based models according to the so-called indirect methods , especially when some proteins are characterized by multiple phosphorylation sites or binding domains, a condition that yields a combinatorial explosion of intermediate chemical complexes and chemical reactions  <cit> .
we designed lassie as a general “black-box” tool able to simulate, in principle, any large-scale reaction-based biochemical system based on mass-action kinetics , given that the available gpu memory is sufficient to accommodate the necessary data structures.
however, considering the difficulty in the manual definition of such massive models, lassie may be adopted as an efficient simulation engine for rule-based modeling tools .
as a matter of fact, rule-based modeling can generate extremely large-scale systems characterized by very long simulation times: lassie may represent an enabling tool to prevent the application of advanced computational investigations of such biological models.
in silico simulations allow to determine the quantitative variation of molecular species amount in time and/or in space, by exploiting either deterministic, stochastic or hybrid algorithms .
in particular, when the concentrations of molecular species is high and the effect of biological noise can be neglected  <cit> , ordinary differential equations  represent the typical modeling approach for cellular processes.
given a model parameterization , the temporal dynamics of the system can be simulated by solving the odes using some numerical integrator, such as euler or runge-kutta methods  <cit> .
unfortunately, odes can be affected by a well-known phenomenon named stiffness  <cit> , which occurs when the system of biochemical reactions is characterized by two well-separated dynamical modes, determined by fast and slow reactions, respectively  <cit> .
stiffness can cause the step-size of integration algorithms to reach extremely small values, thus increasing the overall running time.
to solve this issue, advanced integration methods like lsoda  <cit>  can be exploited, thanks to their capability of efficiently solving stiff systems.
lsoda is able to recognize when a system is stiff and to dynamically select between the most appropriate integration algorithm: the adams methods  <cit>  in the absence of stiffness, and the backward differentiation formulae   <cit>  otherwise.
despite the improvement of efficiency granted by lsoda, the numerical integration of the system of odes can become excessively burdensome when the numbers of reactions and molecular species increase.
lassie overcomes this limitation by distributing over thousands of gpu cores all the calculations required by the numerical integration methods it embeds, therefore paving the way for fast simulations of large-scale and stiff models of cellular processes.
one interesting feature of gpus is that they can have different characteristics, both in terms of resources  and computing power .
kernels’ performances transparently scale on different gpus, since they automatically leverage the additional resources offered by the latest architectures, a characteristic known as transparent scalability.
notably, lassie was designed to be a “black-box” deterministic simulator, not requiring any expertise in mathematical modeling nor any gpu programming skill.
more precisely, given the formalization of a cellular process as a reaction-based model  <cit>  and assuming mass-action kinetics  <cit> , lassie proceeds according to the following workflow:  it automatically generates the system of odes—one ode for each molecular species occurring in the biochemical system—according to the biochemical reactions included in the model;  it automatically derives the jacobian matrix, taking advantage of the symbolic derivation, to apply the bdf;  it executes the numerical integration of the odes by automatically switching between the runge-kutta-fehlberg   <cit>  method in the absence of stiffness and first-order bdf   <cit>  in presence of stiffness.
we point out that lassie is a fully automatic simulator: the user does not need to enter the odes directly.
on the contrary, the input consists in a set of  chemical reactions, specified by means of text files.
the corresponding system of odes is automatically determined according to the mass-action kinetics, making lassie usable without any prior knowledge about odes modeling and integration.
in order to further simplify the execution of simulations, lassie is provided with a user-friendly graphical user interface , whose functioning is described in additional file  <dig>  a comprehensive description of the input files is provided in additional file  <dig>
fig.
<dig> lassie’s graphical user interface that easily allows the user to  open a model,  visualize its set of species, reactions and parameters,  select the output directory,  perform a simulation and  graphically represent the corresponding dynamics
the computational performances of lassie are assessed by measuring the running time required to simulate a set of randomly generated synthetic reaction-based models of increasing size—ranging from  <dig> to  <dig> reactions and species—which is compared to the running time required by a cpu-implementation of lsoda.
moreover, we show the accuracy of lassie by comparing its outcome with lsoda outcome for the simulation of a model of the ras/camp/pka signal transduction pathway in s. cerevisiae  <cit> , which is characterized by stiffness.
we highlight that, in general, the implementation of computational methods able to fully exploit the peculiar architecture of gpus is challenging, since specific programming skills are required and a complete algorithm redesign is often necessary.
for instance, the parallelization on the gpu cores can rely either on a coarse-grained or a fine-grained strategy.
the first strategy allows to simultaneously run a massive number of independent simulations ; on the contrary, the second strategy consists in the parallelization of all the calculations required by a single simulation, an approach that is more suitable for large-scale models.
by virtue of the novel fine-grained parallelization strategy used to implement lassie, our gpu-powered simulator achieves up to 92× speed-up with respect to lsoda.
coarse-grained parallelizations of deterministic simulations were presented in .
the simulators proposed in these works allow to reach a speed-up ranging from 28× to 86× with respect to the corresponding cpu-based simulators.
fine-grained parallelizations of stochastic simulations were presented in  <cit> .
komarov and d’souza proposed gpu-odm  <cit> , a fine-grained simulator of large-scale models based on the stochastic simulation algorithm   <cit> .
this tool uses special data structures and functionalities to efficiently distribute all calculations over the multiple cores of gpus.
these optimizations allow gpu-odm to outperform the most advanced cpu-based implementations of ssa.
komarov et al.
also proposed a gpu-powered fine-grained implementation of τ-leaping  <cit> , an approximate but accurate stochastic algorithm that is, in general, faster than ssa  <cit> .
this tool was shown to be more efficient than its sequential counterpart in the case of extremely large biochemical networks .
notably, to the best of our knowledge, no examples of fine-grained deterministic simulators, such as lassie, have been proposed so far.
lassie was developed using the most widespread gpu computing library, namely, nvidia compute unified device architecture .
cuda allows programmers to exploit the gpus for general-purpose computational tasks .
nevertheless, the direct porting of an application to the gpu is usually unfeasible, so that the full exploitation of the computational power and of the massive parallelism of gpus still represent the main challenges of gpgpu computing.
to exploit the cuda architecture, the programmer implements c/c++ functions , which are loaded from the cpu  to one or more gpus , and replicated in many copies named threads.
cuda organizes threads in three-dimensional structures called blocks, which belong to three-dimensional structures named grids, as shown in fig.
<dig> .
cuda combines the single instruction multiple data  architecture and a flexible multi-threading in order to handle any conditional divergence between threads.
figure  <dig>  also shows a schematic representation of cuda’s memory hierarchy: the global memory , the shared memory , the local memory  and the constant memory .
the global memory is large  but suffers from high access latencies; this problem, anyway, was mitigated thanks to the use of l <dig> cache since the introduction of the fermi architecture.
on the contrary, the constant memory is much smaller  but faster than the global memory, as well as the shared memory ; in particular, the latter should be exploited as much as possible in order to obtain the best performances.
though, the size of the shared memory and its scope restrict the possibility to use it, as only threads belonging to the same block can communicate through the shared memory.
fig.
<dig> threads and memory hierarchy of cuda’s architecture.
left side.
thread organization: the host  launches a single kernel that is executed in multiple threads on the device .
threads  are organized in three-dimensional structures called blocks , which belong to three-dimensional grid .
the programmer must explicitly define the dimensions of blocks and grids.
whenever a kernel is run by the host, the corresponding grid is created by the device which automatically schedules each block on one free streaming multiprocessor available.
this solution allows a transparent scaling of performances on different devices.
moreover, if the machine is equipped with more than one gpu, it is also possible to distribute the workload by launching the kernel on each gpu.
right side.
memory hierarchy: in cuda there are many different memories with different scopes.
each thread has two different kind of private memory: registers and local memories.
threads belonging to the same block can communicate through the shared memory, which has low access latency.
the global memory suffers from high access latencies but it is accessible to all threads and it is cached since the introduction of the fermi architecture.
also the texture and the constant memory are equipped with a cache as well, and all threads can read from these two memories
given the peculiar features of cuda architecture, gpu programmers should be able to optimize both threads partitioning and memory usage, as well as to redesign the algorithm with appropriate kernels, in order to fully leverage the computational power of these multi-core devices.
for instance, in the implementation of lassie, the shared memory is not used because several blocks are exploited to solve the system of odes, and threads do not communicate data with each other.
moreover, the data structures employed by lassie are larger than the total size of the shared memory, thus preventing the possibility to exploit it.
in what follows, we show how gpu programming and cuda features—including built-in support for vector types, which extend the standard c data types to vector—have been exploited to optimize the execution workflow of lassie.
the paper is structured as follows.
in the next section we briefly introduce the formalism of reaction-based models and provide a general description of lassie’s implementation.
then, we discuss the computational performance of lassie, showing the speed-up it achieves with respect to lsoda for the simulation of reaction-based models of different sizes.
we also analyze how the number of reactions and the number of species affect the performances of lassie.
we conclude the work with some final remarks about cuda’s architecture and lassie, proposing future improvements of the simulator.
lassie is available on the github repository https://github.com/aresio/lassie.
methods
lassie is designed to be a “black-box” deterministic simulator, created to be easily used without any gpu programming or odes modeling skills.
in this section we describe how lassie allows to perform deterministic simulations of large-scale biochemical models, distributing all required calculations on the cores of the gpu.
it is worth noting that the parallelization strategy exploited by lassie represents one of the novelties of this work, and allowed to achieve the remarkable performance results presented in the next sections.
in particular, lassie has been developed to solve systems of coupled odes specified in the form dxdt=f, where x≡x represents the vector of concentration values at time t of all chemical species occurring in the system.
reaction-based models and odes generation
reaction-based modeling is a mechanistic, quantitative and parametric formalism to describe and simulate networks of biochemical reactions  <cit> , which was exploited to analyze different signal transduction pathways .
a reaction-based model is defined by specifying the set of n molecular species {s
<dig> …,s
n} and the set of m biochemical reactions {r
<dig> …,r
m} which appear in the cellular process under investigation  <cit> .
a generic reaction is described as follows:
<dig> ri:∑j=1naijsj→ki∑j=1nbijsj,i= <dig> …,m,
where a
ij, bij∈ℕ are the stoichiometric coefficients and ki∈ℝ+ is the kinetic constant associated with r
i.
the set of reactions {r
<dig> …,r
m} can be written compactly in the matrix-vector form as→kbs, where s=[s
1⋯s
n]t is the n-dimensional column vector of molecular species, k=[k
1⋯k
m]t is the m-dimensional column vector of kinetic constants, and a,b∈n
m×n are the so-called stoichiometric matrices whose  elements i,j and i,j correspond to the stoichiometric coefficients a
ij and b
ij of the reactants and the products of all reactions, respectively.
since a reaction simultaneously involving more than two reactants has a probability to take place almost equal to zero, here we consider only first and second-order reactions .
for this reason, the matrices a and b are sparse.
given an arbitrary reaction-based model and assuming the law of mass-action  <cit> , it is possible to derive the corresponding system of coupled odes that describes the variation in time of the species concentrations.
specifically, by denoting the concentration of species s
j at time t as x
j, where x
j∈r
≥ <dig> for j= <dig> …,n, the system of coupled odes can be obtained as follows:
<dig> dxdt=t,
where x is the n-dimensional vector of concentration values at time t , the symbol ∘ denotes the entry-by-entry matrix multiplication, and x
a denotes the vector-matrix exponentiation form  <cit> .
formally, x
a is a m-dimensional vector whose i-th component is given by x1ai1⋯xnain, for i= <dig> …,m.
we highlight that each ode appearing in eq.
<dig> is a polynomial function, consisting in at least one monomial that is associated with a specific kinetic constant.
data structures and cuda memory usage
given a reaction-based model as input, lassie automatically generates the systems of odes according to eq.
<dig> and encodes the matrices a and h=t as two arrays of short <dig> cuda vector types, named v
a and v
h, respectively.
cuda vector types are multi-dimensional data ranging from  <dig> to  <dig> components, addressed by .x, .y, .z, and .w.
since the matrices a and h are sparse, lassie uses compressed data structures created by removing all zero elements from a and h, in order to save memory and avoid unnecessary readings from the global memory.
namely, let h
ji be the element of h at row j and column i, and a
ij the element of a at row i and column j, for i= <dig> …,m and j= <dig> …,n.
for each non-zero element of h, we store into the .x and .y components of v
h the values j and i, respectively; the .z component of v
h is used to store the element h
ji, while the .w component stores the index of the kinetic constant associated with that monomial.
similarly, for each non-zero element of a, the .x and .y components of v
a contain the values i and j, respectively.
the value a
ij is stored into the .z component of v
a, while the .w component is left unused.
note that we exploited the short <dig> cuda vector type rather than the short <dig> cuda vector type, because the former is 8-aligned and requires a single instruction to fetch a whole entry, while the latter is 2-aligned and thus takes three memory operations to read each entry.
in order to parse these arrays inside the gpu, we use two additional arrays of short <dig> cuda vector types, named o
h and o
a, which store the offsets used to correctly read the entries of the v
h and v
a structures, respectively.
the .x and .y components of each row of o
h contain, respectively, the first index and the last index to access the v
h structure.
each thread uses its own pair of indexes to read the rows of the v
h structure between the first index and the last one.
similarly, o
a stores the indexes that allow to correctly access the v
a structure.
finally, the values of the kinetic constants are stored into an array of type double, named k. figure  <dig> shows an example of the matrix encoding used in lassie.
fig.
<dig> example of matrix encoding to automatically generate an ode using lassie.
all terms of the polynomial function describing the ode of species x
<dig> given at the top of the figure are encoded in the components of the data structures o
h, v
h, o
a, v
a and k, as detailed hereby.
notice that only the data structures components with solid borders are used to automatically generate the ode; the various terms appearing in the ode are represented with corresponding colors in the data structure components.
matrix encoding starts from matrix o
h. each thread j, for j= <dig> …,n− <dig>  reads the values stored in the .x and .y components of o
h .
in this example, we consider species x
<dig> that corresponds to thread  <dig>  each thread fetches the values in v
h, starting from the row indicated by the value stored in the .x component of o
h, up to the row corresponding to the value stored in the .y component.
in this example, thread  <dig> in matrix o
h reads the values contained in the first two rows—i.e., rows  <dig> and 1—in matrix v
h. each row of v
h encodes a monomial of an ode: the .x component is not used; the .y components  indicate the row numbers of the o
a structure that each thread must read; the .z components  indicate the sign and the coefficient of the monomial; the .w components  indicate the positions of the array k containing the values of the kinetic constants corresponding to the reactions that the threads are parsing.
in this example, the .z and .w components of v
h allows to derive the coefficients −1k
<dig> and +1k
<dig> for the first and the second term of the ode, respectively.
afterwards, as in the case of o
h, each thread fetches the values in v
a, starting from the row indicated by the value stored in the .x component of o
a, up to the row corresponding to the value stored in the .y component of o
a. the values stored in the .y  and .z  components of v
a correspond to the indexes of the species and the stoichiometric coefficients, respectively, while the .x and .w components of v
a are left unused.
in this example, row  <dig> in matrix o
a reads the values stored in rows  <dig> and  <dig>  of matrix v
a, generating the factors  <dig> in the first term of the ode, while row  <dig> in matrix o
a reads the values stored in row  <dig>  of matrix v
a, generating the factor  <dig> in the second term of the ode.
therefore, in this example, the matrix encoding overall generates the ode of species x
<dig> consisting in the sum of two polynomial terms: −k
111+k
21
thanks to these cuda structures, we obtain a twofold performance improvement:  at the instruction level, a single instruction is enough to either load or store a multi-word vector.
so doing, the total instruction latency for a particular memory transaction is lower and also the bytes per instruction ratio is higher;  at the memory controller level, by using vector types a transfer request from a warp has a larger net memory throughput per transaction, yielding a higher bytes per transaction ratio.
with a fewer number of transfer requests, the memory controller is able to reduce contentions producing a higher overall memory bandwidth utilization.
the only limitation due to short data type is that indices are limited to 22×8− <dig>  which means that lassie cannot simulate systems larger than  <dig>  <dig> chemical species and reactions.
execution workflow and cuda kernels
once that the system of odes is generated by reading the input files  and appropriately stored according to the cuda vector types, lassie solves it by automatically switching between the runge-kutta-fehlberg  method  <cit>  in the absence of stiffness, and the backward differentiation formulae  methods  <cit>  in presence of stiffness.
the integration of the systems of odes is carried out from an initial time instant t
<dig>  up to a given maximum simulation time t
max.
in order to reproduce the dynamics of the cellular process described by the odes, the concentration values of the molecular species appearing in the reaction-based model are saved at specified time steps within the interval [t
<dig> t
max] .
lassie’s workflow consists in  <dig> distinct phases, as represented in fig.
<dig>  note that phases p
<dig>  p
<dig> and p
<dig> are executed by the host , while p
<dig>  p
<dig> and p
<dig> are executed by the device .
overall, phases p
<dig>  p
<dig> and p
<dig> rely on  <dig> different lightweight kernels, which were specifically developed to fully leverage the parallel architecture of the gpu for the implementation of the aforementioned numerical integration methods.
we describe hereafter the main design and implementation choices of each phase and their related cuda kernels, which result in a novel parallelization strategy with respect to state-of-the-art methodologies .
fig.
<dig> simplified scheme of lassie workflow.
the data structures used to encode the system of odes are generated in phase p
<dig>  in phase p
<dig>  if the current simulation time t corresponds to a specified sampling time instant, then the current concentration values of all molecular species are saved; otherwise, the execution proceeds to the next phase.
in phase p
<dig> each thread derives and solves the corresponding ode by exploiting the rkf method, while in phase p
<dig> the rkf solutions are verified:  if the rkf solutions are rejected, then the integration step-size dt is reduced and phase p
<dig> is executed again;  if rkf solutions are rejected but the integration step-size dt is too small, then phase p
<dig> is executed and the system of odes is solved using the bdf methods;  if the rkf solutions are accepted, the termination criterion is verified during phase p
<dig>
phase
p
1
.
it implements the generation of all data structures used to encode the odes, as described in the previous section.
this phase is executed on the host.
phase
p
2
.
it is used to sample and save the system dynamics and it is implemented by means of a single cuda kernel .
in particular, if the current simulation time t corresponds to one of the specified sampling time instants, lassie saves the concentration values of  all molecular species into an array defined on the gpu.
otherwise, the execution proceeds to the next phase.
phase
p
3
.
it implements the rkf method  <cit> , an explicit integration algorithm with variable step-size used by each thread j to solve the j-th ode, for j= <dig> …,n− <dig>  this phase is implemented as  <dig> cuda kernels.
during this phase, two different approximated states u and w of the state x of the system are generated at each step, thanks to the evaluation of six supplementary values l
<dig> …,l
<dig> .
to evaluate the accuracy of u and w at the current step-size dt, lassie exploits a user-defined vector tolerance ε∈r
n , and two additional arrays, e
r,δ∈r
n, defined as follows:
<dig> er=|w−u|dt,δ= <dig> εer <dig>
if e
r
j≤ε
j for all j= <dig> …,n, then u is accepted as new state of the system, that is, x=u; otherwise, the solutions u and w are rejected and recalculated by using a new step-size.
the new step-size is computed as d
t=d
t· min{δ
<dig> …,δ
n}, being δ
<dig> …,δ
n the components of vector δ .
overall, phase p
<dig> is implemented by means of the following kernels:
kernel
k
2: used to evaluate each ode at the current state x of the system;
kernels
k
<dig> – k
8: each thread j, for j= <dig> …,n− <dig>  computes the components l
1j,…,l
6j of l
<dig> …,l
<dig>  by invoking kernel
k
2;
kernel
k
9: each thread j, for j= <dig> …,n− <dig>  computes the components w
j and u
j of the approximated states u and w, respectively;
kernel
k
10: each thread j, for j= <dig> …,n− <dig>  calculates the components e
r
j and δ
j of e
r and δ, respectively.
phase
p
4
.
it is used to verify the rkf solutions calculated during phase p
<dig> and, accordingly, to choose the next phase to be executed:  if the solutions are rejected and the new step-size dt is acceptable , phase p
<dig> is executed again exploiting a smaller step-size dt;  if the solutions are rejected and the new step-size dt becomes too small , lassie executes phase p
5;  if all solutions do not violate the specified rkf-tolerance vector ε, then lassie executes phase p
<dig>
note that point  implicitly states that the system of odes is considered to be stiff, so that lassie automatically switches to phase p
<dig>  where bdf methods are used for the numerical integration.
phase p
<dig> is executed on the host.
phase
p
5
.
it implements the bdf methods, the most widely used implicit multi-step numerical integration algorithms  <cit> .
lassie switches to this phase if and only if the rkf solutions u and w evaluated during phase p
<dig> are rejected, and the rkf step-size dt becomes smaller than ε
s.
the general formula for a bdf can be written as
<dig> ∑i=0qαix=dtβ0f),
where the coefficients α
i  and β
<dig> are chosen according to the order q of bdf  <cit> , and dt is user-defined.
note that, for q> <dig>  the absolute stability region of the resulting bdf methods is too small and such bdfs are numerical unstable  <cit> .
therefore, bdfs with an order q greater than  <dig> are not used.
since each bdf is an implicit method, at each time step it requires the solution of a nonlinear system of equations, which can be solved by using the iterative newton—raphson method  <cit> .
this algorithm allows to find successively better approximations z of the zeros of a real-valued function f= <dig> by using the derivative of f, and it is repeated until a sufficiently accurate value is reached.
this idea can be extended to a system of nonlinear equations, by using the jacobian matrix j) of f), which is the matrix of all first-order partial derivatives.
since the evaluation of the jacobian matrix at each iteration is computationally expensive, lassie actually exploits:  a modified newton—raphson method  <cit> ;  the lu factorization method  <cit>  .
during phase p
<dig>  the newton-raphson method is iterated until a user-defined maximum number of iterations m
a
x
it is reached, or a sufficiently accurate value is achieved .
overall, phase p
<dig> is implemented by means of the following kernels:
kernel
k
11: each thread j, for j= <dig> …,n− <dig>  derives the j-th row of the jacobian matrix and evaluates it on the current state of the system x;
kernel
k
12: the jacobian matrix is transposed in order to exploit the lu factorization method ;
kernels
k
<dig> – k
18: based on the order q of the bdf, lassie invokes one of these kernels  to calculate the known terms of the linear system;
kernels
k
<dig> – k
24: each kernel
k
, q= <dig> …, <dig>  performs the calculations of the q-th order bdf;
kernel
k
25: it updates the iteration vector needed to execute the newton-raphson method .
phase
p
6
.
it is used to verify the termination criterion: if the maximum time t
max is reached, then the simulation ends.
on the contrary, the execution iterates from phase p
<dig>  this phase is executed on the host.
all the temporary results computed by lassie are stored on the gpu, since data transfers between the host and the device are very time consuming.
for the same reason, the output data  are transferred to the host as soon as the whole simulation is completed.
RESULTS
in this section we compare the computational performance of lassie against lsoda  <cit> , which is generally considered one of the best numerical integration algorithms for deterministic simulations of biological systems, thanks to its capability of dealing with stiff and non-stiff systems.
in particular, we exploited the lsoda implementation provided by scipy library  <cit>  , written in c language.
lassie was run on a machine with a gpu nvidia geforce titan gtx, based on the kepler architecture and equipped with 2× <dig> streaming multiprocessors for a total of  <dig> cores  and a theoretical peak processing power of  <dig>  tflops in double precision.
instead, lsoda was run on galileo, a supercomputer created by the italian consortium cineca.
galileo consists of  <dig> compute nodes, each one equipped with  <dig> cpus octa-core intel xeon haswell e5- <dig> v <dig>  for a total of  <dig> cores, and  <dig> gb of ram.
each cpu is capable of about  <dig> gflops in double precision.
in our tests, we exploited one node with  <dig> gb of ram distributed over  <dig> cores.
the computational performance was evaluated by simulating a set of synthetic reaction-based models of increasing size, that is, having a number of reactions and species m×n arbitrarily chosen in the range from 64× <dig> to 8192× <dig>  the models were generated considering the methodology used in  <cit> , which was modified in order to randomly sample the initial concentration of each species with a uniform distribution in the range [ <dig> ), and the kinetic constant of each reaction with a logarithmic distribution in the range [10− <dig> ).
for each model size m×n, we generated and simulated  <dig> different synthetic reaction-based models to the aim of measuring the average running time of both lassie and lsoda.
the simulation of each reaction-based model was performed multiple times, using different settings for the sampling of the time-series.
specifically, in each repetition, we saved either  <dig> , <dig>  or  <dig> samples of the system dynamics of all chemical species, at regular intervals.
all simulations were halted at time t
max= <dig> .
all simulations were executed—independently from the size of the model and the number of samples saved—by setting the following parameters of lassie:
tolerance of rkf method ε
j = 10− <dig>  j= <dig> …,n;
first–order bdf method ;
bdf integration step d
t= <dig> ;
tolerance of newton-raphson method ε
nr = 10−6;
maximum number of iterations allowed during each call of the newton-raphson method m
a
x
it=104;
initial integration step of rkf method equal to 10−3;
tolerance value to switch between rkf and backward euler methods ε
s=10− <dig>
the following parameters of lsoda were used to run the simulations:
relative tolerance equal to 10−6;
absolute tolerance equal to 10−12;
maximum number of internal steps equal to  <dig>
fig.
<dig> speed-up values  achieved by lassie with respect to lsoda for the simulation of synthetic models of increasing size, having a number of reactions and of species m×n  and characterized by an increasing number of sampling time instants of the system dynamics .
when the value of the speed-up is greater than one, lassie is faster than lsoda and vice versa
m×n
lsoda
lassie
speed-up
lsoda
lassie
speed-up
lsoda
lassie
speed-up
lsoda
lassie
speed-up
lsoda
lassie
speed-up
na
na
na
na
na
na
na
na
na
na
*maximum speed-up value
fig.
<dig> comparison between the average running time required by lassie  and lsoda  to simulate  <dig> instances of models characterized by  <dig> reactions and  <dig> species ,  <dig> reactions and  <dig> species ,  <dig> reactions and  <dig> species , saving different numbers of sampling time instants of the dynamics.
note that the y-axes are in logarithmic scale
as an additional test, we investigated whether the relationship between the number of reactions and the number of species could affect the overall performances of lassie.
as the number of chemical species corresponds to the number of odes, the length of each ode is roughly proportional to the number of reactions.
since gpus have a lower clock frequency than cpus , each gpu core is slower than the cpu core to perform a single instruction <dig>  for this reason, in order to obtain the highest performances, the calculations on the gpu should be spread across threads as much as possible, while the number of operations performed by each thread should be reduced.
indeed, as reported in table  <dig> and shown in fig.
<dig>  when the number of chemical species involved in a model is greater than the number of reactions, lassie achieves better performances than those obtained in the case of models with a number of chemical species smaller than the number of reactions.
for instance, considering the models with m×n equal to 171× <dig>  the running time of lassie is smaller than in the case of the models with size 512× <dig>  irrespective of the number of samples of the system dynamics, thanks to the higher number of threads that are concurrently launched on the gpu in the first case.
this is in general valid in all cases with the exception of the models characterized by  <dig> chemical species with  <dig> and  <dig> samples of the system dynamics.
here, the average running time of lassie is greater than in the case of models with  <dig> reactions, since the required number of accesses to the high-latency global memory of the gpu impairs the performances of the simulations.
fig.
<dig> running time  of lassie for the simulation of synthetic models of increasing size, having a number of reactions and of species m×n , with m≠n, and characterized by an increasing number of sampling time instants of the system dynamics
m×n
lassie
lassie
lassie
lassie
lassie
in order to assess the scalability of lassie, and of cuda applications in general, we executed additional tests on different gpus.
figure  <dig> shows a comparison of lassie’s performance using three different gpu models : a notebook video card , the nvidia geforce titan z used throughout the paper , and a tesla-class gpu .
to compare the speed-up provided by these gpus we generated  <dig> different synthetic models  and calculated the average running time.
fig.
<dig> comparison of the average running times for the simulation of  <dig> synthetic models characterized by three different sizes, executed with different gpus: a notebook gpu nvidia geforce 960m ; a nvidia geforce titan z ; a tesla-class gpu nvidia k20c
our results highlight the importance of two distinct factors on lassie’s performances: the gpu’s clock frequency and the amount of available resources .
as a matter of fact, despite the lower amount of cuda cores, the geforce 960m turns out to be competitive on models of moderately large size thanks to its higher clock rate, with respect to the titan z and the k20c.
when the odes largely outnumber the available cores , the geforce 960m is no longer competitive.
this is an example of transparent scalability of cuda applications: the threads are automatically distributed over the available cores, improving the overall performances, without any user intervention.
moreover, as described in the background section, threads are organized in blocks that are scheduled on the available multi-processors.
thanks to this characteristic, when the overall number of threads outnumbers the available cores, cuda automatically creates a queue of blocks that are scheduled on the streaming multi-processors as soon as they become available for computation.
thus, lassie can, in principle, simulate any model on any gpu, as long as there is enough memory to store the data structures.
the tesla k20c is characterized by a large amount of cores that, in the case of 1024× <dig> models, are fully exploited only during the simulation of the stiff parts of the dynamics.
for the remaining parts of the simulation, half of its cores are actually used for computation with a slower clock rate with respect to the clock rate of the geforce gpus.
moreover, tesla cards exploit error correcting codes  on memories, ensuring additional checks of correctness to the data against potential corruption from electrical or magnetic interference, at the price of a significant overhead  <cit> .
the ecc was enabled during all tests, partly explaining the reduced performance of the tesla k20c on very large-scale models with respect to the titan z.
we assessed the accuracy of lassie by simulating the dynamics of the model of the ras/camp/pka signaling pathway in yeast presented in  <cit> , and comparing the outcome of lassie with the result of the simulation performed with lsoda.
we also investigated the influence of lassie parameters  on the running times and quality of the simulated solutions, by exploiting a model representing a chain of isomerizations.
the accuracy results—which show an identical dynamics with respect to lsoda using default settings—are presented in the additional file  <dig>
as a final remark, we highlight that a fair comparison of gpus and cpus is a difficult task, in general, due to their deep architectural differences.
the theoretical peak performances of both architectures are difficult to achieve: indeed, developers must implement code to the aim of maximizing the parallelism and the occupancy of the multi-processors, adhering as much as possible to the underlying simd computational model in the case of the gpu and exploiting vector instructions in the case of the cpu.
however, gpus allow the temporary divergence of the execution flow of threads, that is, a part of the threads can execute different portions of the code .
when this situation occurs, some threads get stalled waiting for reconvergence.
this mechanism provides the programmer with a certain degree of freedom to abandon the simd paradigm, but at the same time it can potentially lead to the complete serialization of the execution affecting the overall performances.
hence, conditional branches should be avoided as much as possible.
we also highlight that the usage of registers and shared memory influences the occupancy of the gpu, as these resources are scarce on each streaming multiprocessor.
all these circumstances can prevent the achievement of the peak computational power of a gpu.
to this aim, we developed kernels that maximize the parallelism and the occupancy of the multi-processors avoiding threads divergence as much as possible.
moreover, we optimized data structures to store the matrices a and h that encode the system of odes, and cuda vector types that allow to increase the memory throughput and to reduce the number of memory accesses, all precautions that explain the performance boost achieved with lassie.
CONCLUSIONS
in this work we presented lassie, a gpu-powered simulator of large-scale biochemical systems based on mass-action kinetics.
lassie is a “black-box” simulator able to automatically convert reaction-based models of biological systems into the corresponding systems of odes.
reaction-based models defined according to the law of mass-action do not hinge upon the use of any approximate kinetics functions , which are frequently used in systems biology for the definition of mathematical models based on differential equations.
although michaelis-menten kinetics or hill functions can be useful in biological modeling, they rely on chemical assumptions that are valid only in certain conditions  <cit> .
therefore, the reason why we rely on mass-action based models is manifold.
on the one hand, since the biological function and biochemical kinetics of all molecular species and all reactions appearing in the model are not approximated nor lumped together in any way, they can be analyzed independently from each other.
as a consequence, this allows to determine the influence of every single species and reaction on the overall functioning of the system.
on the other hand, the law of mass-action allows to derive a first order ode for each species appearing in the model: it is worth noting that such ode is a polynomial function that describes how the concentration of that species changes in time, according to all the reactions where it appears either as reactant or product  <cit> .
the presence of polynomial functions simplifies the symbolic derivation that is needed to calculate the jacobian matrix associated with the odes and exploited by the bdf.
in addition, as described in the methods section, polynomials can be efficiently encoded in the memory and parsed gpu-side.
as a result, all gpu threads can perform the same task , strongly reducing warps’ divergence and the consequent stalling of threads due to serialization, a circumstance that would instead happen if each thread calculated an ode characterized by an arbitrary kinetics.
in order to solve systems of odes characterized by stiffness, lassie automatically switches between the rkf and the bdf integration methods.
lassie’s execution flow is partitioned into  <dig> cuda kernels, overall distributing the calculations over the available cores in order to fully exploit the massive parallel capabilities of modern gpus, therefore achieving a relevant reduction of the running time in case of large-scale models.
in order to assess the computational performance of lassie, we performed a set of simulation tests using synthetic reaction-based models of increasing size, and we compared lassie’s running time with respect to the lsoda numerical integration algorithm implemented in the scipy library.
the break-even between the performances of lassie and lsoda was observed when both the numbers of reactions and chemical species is in between  <dig> and  <dig>  this result indicates that, for biological systems consisting in more than  <dig> reactions and  <dig> species, the gpu-powered simulator becomes more convenient than the lsoda algorithm running on cpu.
indeed, in the case of large-scale models, characterized by  <dig> reactions and  <dig> species, we obtained a considerable 92× speed-up.
moreover, thanks to its smaller memory footprint with respect to lsoda, lassie allows the simulation of even larger models, taking just an average of  <dig>  s to simulate models characterized by  <dig> reactions and  <dig> species.
on the contrary, lsoda did not allow the simulation of models of this size on the computer we used for the tests, as it crashed because of its very high memory footprint.
we also highlight that copasi  <cit> , one of the most used software in systems biology, requires in general longer execution times with respect to the scipy implementation of lsoda exploited in this work.
in addition, copasi fails when trying to simulate large-scale models.
we provide an example of such model—characterized by  <dig> reactions and  <dig> species—as sbml file in the github repository.
bdfs are the most used integration algorithms to solve systems of odes in case of stiffness.
the first–order bdf is a single-step implicit integration method, meaning that the next state of the system depends only on the current state of the system.
higher–order bdfs are multi-step methods, so that the next state of the system relies on multiple previous states of the system .
this implies that the integration step-size should be the same for all previous states to ensure the correctness of the solution.
for this reason, lsoda uses the multi-step adams methods as explicit methods in addition to bdfs.
conversely, lassie uses the rkf method, which is a single-step explicit algorithm with variable step-size, and the backward euler method.
other single-step implicit methods belonging to the family of runge-kutta methods exist  <cit> , the most known being the families of lobatto and radau methods  <cit> .
these methods have been proven to be suitable for stiff systems, thanks to their accuracy and stability  <cit> .
as a future development of lassie, we will investigate the feasibility and efficacy of replacing the backward euler and, more generally, the bdfs with implicit runge-kutta methods  <cit> , in particular lobatto and radau methods.
in order to fully exploit the cuda architecture, the memory hierarchy must be exploited as much as possible.
because of the peculiar sequential structure of both explicit and implicit integration algorithms, lassie’s kernel are lightweight and rarely reuse any variables.
for this reason, the current implementation only leverages the global memory  and registers to manipulate the mutable data.
the shared memory has not been exploited in any way, leaving room for potential future improvements of performances.
however, on the gpus where the l <dig> cache and the shared memory share the same resources, cuda allows to express a preference to assign a larger amount of memory to the caching mechanisms.
this functionality is enabled by default in lassie using the cuda cudafuncsetcacheconfig primitive, executed with the cudafunccachepreferl <dig> argument.
also the constant memory has not been used, since all data structures are larger than the total size of this memory.
in a future release, we plan to leverage these memories, for instance to store the array of the kinetic constants and the structures containing the offsets used to correctly decode the odes.
lassie currently exploits only a single gpu, even on multi-gpu systems; as a future improvement of this work we plan to extend it in order to support multi-gpu systems, to further increase the size of the models to be simulated.
finally, an additional goal in the development of lassie is to integrate and accelerate the investigation of rule-based models.
in the next future, we plan to develop a set of tools that will leverage the functionalities offered by rule-based modeling frameworks , to convert a rule-based model into a set of reactions.
although rule-based tools already provide internal simulation methods , lassie can represent a valuable alternative for large-scale models, enabling the investigation of more detailed biological systems, paving the way to potential new discoveries in systems biology.
lassie will be also integrated in cosys, a free web-based platform for systems biology investigation available at http://www.sysbio.it/cosys  <cit> .
endnote
<dig> the advances in gpu’s technology will progressively reduce this gap.
in middle 2016—with the introduction of the novel pascal architecture and the  <dig> nm finfet manufacturing process—nvidia presented a gpu with a clock frequency of  <dig>  ghz that, theoretically, is expected to double lassie’s performances.
additional files
additional file  <dig> lassie graphical user interface.
additional file  <dig> lassie input files and command line arguments.
additional file  <dig> implementation of cuda kernels for lassie execution workflow.
additional file  <dig> simulation accuracy of lassie.
abbreviations
bdfbackward differentiation formulae
cpucentral processing unit
cudacompute unified device architecture
eccerror correcting codes
gpgpugeneral-purpose gpu
gpugraphics processing unit
lassielarge-scale simulator
odeordinary differential equation
rkfrunge-kutta-fehlberg
simdsingle instruction multiple data
ssastochastic simulation algorithm

 systems biology is a multidisciplinary research field relying on the cross-talk between mathematical, computational and experimental tools to investigate the functioning of complex biological systems, and to predict how they might behave in both physiological and perturbed conditions. to this aim, different computational methods—e.g., parameter estimation, sensitivity analysis or reverse engineering —are usually exploited to define or calibrate the mathematical model that describes the system of interest. these methods require the execution of a large number of simulations, each one generally corresponding to a distinct model structure or parameterization, that is, to a different set of molecular interactions or to different initializations of the species amounts and/or reaction constants. as a result, the computational burden required by these computational analyses can rapidly overtake the capabilities of central processing units, therefore limiting in-depth computational investigations to small-scale models consisting in a few tens of reactions and molecular species at most. general-purpose graphics processing units can be exploited to overcome these drawbacks. indeed, they are parallel multi-core co-processors that are drawing an ever-growing attention by the scientific community, since they give access to tera-scale performances on common workstations. as such, they can markedly decrease the running times required by traditional cpu-based software, still maintaining low-costs and energetic efficiency. as a matter of fact, in the latter years gpus have been widely adopted as an alternative approach to classic parallel architectures for the parallelization of computational methods in systems biology, computational biology and bioinformatics. in this work we propose lassie, a novel gpu-accelerated software designed to simulate large-scale reaction-based models of cellular processes, consisting in hundreds or thousands of reactions and molecular species. an example of killer-application of lassie would consist in the simulation of rule-based models according to the so-called indirect methods, especially when some proteins are characterized by multiple phosphorylation sites or binding domains, a condition that yields a combinatorial explosion of intermediate chemical complexes and chemical reactions. we designed lassie as a general “black-box” tool able to simulate, in principle, any large-scale reaction-based biochemical system based on mass-action kinetics, given that the available gpu memory is sufficient to accommodate the necessary data structures. however, considering the difficulty in the manual definition of such massive models, lassie may be adopted as an efficient simulation engine for rule-based modeling tools. as a matter of fact, rule-based modeling can generate extremely large-scale systems characterized by very long simulation times: lassie may represent an enabling tool to prevent the application of advanced computational investigations of such biological models. in silico simulations allow to determine the quantitative variation of molecular species amount in time and/or in space, by exploiting either deterministic, stochastic or hybrid algorithms. in particular, when the concentrations of molecular species is high and the effect of biological noise can be neglected, ordinary differential equations represent the typical modeling approach for cellular processes. given a model parameterization, the temporal dynamics of the system can be simulated by solving the odes using some numerical integrator, such as euler or runge-kutta methods. unfortunately, odes can be affected by a well-known phenomenon named stiffness, which occurs when the system of biochemical reactions is characterized by two well-separated dynamical modes, determined by fast and slow reactions, respectively. stiffness can cause the step-size of integration algorithms to reach extremely small values, thus increasing the overall running time. to solve this issue, advanced integration methods like lsoda can be exploited, thanks to their capability of efficiently solving stiff systems. lsoda is able to recognize when a system is stiff and to dynamically select between the most appropriate integration algorithm: the adams methods in the absence of stiffness, and the backward differentiation formulae otherwise. despite the improvement of efficiency granted by lsoda, the numerical integration of the system of odes can become excessively burdensome when the numbers of reactions and molecular species increase. lassie overcomes this limitation by distributing over thousands of gpu cores all the calculations required by the numerical integration methods it embeds, therefore paving the way for fast simulations of large-scale and stiff models of cellular processes. one interesting feature of gpus is that they can have different characteristics, both in terms of resources and computing power. kernels’ performances transparently scale on different gpus, since they automatically leverage the additional resources offered by the latest architectures, a characteristic known as transparent scalability. notably, lassie was designed to be a “black-box” deterministic simulator, not requiring any expertise in mathematical modeling nor any gpu programming skill. more precisely, given the formalization of a cellular process as a reaction-based model and assuming mass-action kinetics, lassie proceeds according to the following workflow: it automatically generates the system of odes—one ode for each molecular species occurring in the biochemical system—according to the biochemical reactions included in the model; it automatically derives the jacobian matrix, taking advantage of the symbolic derivation, to apply the bdf; it executes the numerical integration of the odes by automatically switching between the runge-kutta-fehlberg method in the absence of stiffness and first-order bdf in presence of stiffness. we point out that lassie is a fully automatic simulator: the user does not need to enter the odes directly. on the contrary, the input consists in a set of chemical reactions, specified by means of text files. the corresponding system of odes is automatically determined according to the mass-action kinetics, making lassie usable without any prior knowledge about odes modeling and integration. in order to further simplify the execution of simulations, lassie is provided with a user-friendly graphical user interface, whose functioning is described in additional file a comprehensive description of the input files is provided in additional file fig. lassie’s graphical user interface that easily allows the user to open a model, visualize its set of species, reactions and parameters, select the output directory, perform a simulation and graphically represent the corresponding dynamics the computational performances of lassie are assessed by measuring the running time required to simulate a set of randomly generated synthetic reaction-based models of increasing size—ranging from to reactions and species—which is compared to the running time required by a cpu-implementation of lsoda. moreover, we show the accuracy of lassie by comparing its outcome with lsoda outcome for the simulation of a model of the ras/camp/pka signal transduction pathway in s. cerevisiae, which is characterized by stiffness. we highlight that, in general, the implementation of computational methods able to fully exploit the peculiar architecture of gpus is challenging, since specific programming skills are required and a complete algorithm redesign is often necessary. for instance, the parallelization on the gpu cores can rely either on a coarse-grained or a fine-grained strategy. the first strategy allows to simultaneously run a massive number of independent simulations; on the contrary, the second strategy consists in the parallelization of all the calculations required by a single simulation, an approach that is more suitable for large-scale models. by virtue of the novel fine-grained parallelization strategy used to implement lassie, our gpu-powered simulator achieves up to 92× speed-up with respect to lsoda. coarse-grained parallelizations of deterministic simulations were presented in. the simulators proposed in these works allow to reach a speed-up ranging from 28× to 86× with respect to the corresponding cpu-based simulators. fine-grained parallelizations of stochastic simulations were presented in. komarov and d’souza proposed gpu-odm, a fine-grained simulator of large-scale models based on the stochastic simulation algorithm. this tool uses special data structures and functionalities to efficiently distribute all calculations over the multiple cores of gpus. these optimizations allow gpu-odm to outperform the most advanced cpu-based implementations of ssa. komarov et al. also proposed a gpu-powered fine-grained implementation of τ-leaping, an approximate but accurate stochastic algorithm that is, in general, faster than ssa. this tool was shown to be more efficient than its sequential counterpart in the case of extremely large biochemical networks. notably, to the best of our knowledge, no examples of fine-grained deterministic simulators, such as lassie, have been proposed so far. lassie was developed using the most widespread gpu computing library, namely, nvidia compute unified device architecture. cuda allows programmers to exploit the gpus for general-purpose computational tasks. nevertheless, the direct porting of an application to the gpu is usually unfeasible, so that the full exploitation of the computational power and of the massive parallelism of gpus still represent the main challenges of gpgpu computing. to exploit the cuda architecture, the programmer implements c/c++ functions, which are loaded from the cpu to one or more gpus, and replicated in many copies named threads. cuda organizes threads in three-dimensional structures called blocks, which belong to three-dimensional structures named grids, as shown in fig.. cuda combines the single instruction multiple data architecture and a flexible multi-threading in order to handle any conditional divergence between threads. figure also shows a schematic representation of cuda’s memory hierarchy: the global memory, the shared memory, the local memory and the constant memory. the global memory is large but suffers from high access latencies; this problem, anyway, was mitigated thanks to the use of l cache since the introduction of the fermi architecture. on the contrary, the constant memory is much smaller but faster than the global memory, as well as the shared memory; in particular, the latter should be exploited as much as possible in order to obtain the best performances. though, the size of the shared memory and its scope restrict the possibility to use it, as only threads belonging to the same block can communicate through the shared memory. fig. threads and memory hierarchy of cuda’s architecture. left side. thread organization: the host launches a single kernel that is executed in multiple threads on the device. threads are organized in three-dimensional structures called blocks, which belong to three-dimensional grid. the programmer must explicitly define the dimensions of blocks and grids. whenever a kernel is run by the host, the corresponding grid is created by the device which automatically schedules each block on one free streaming multiprocessor available. this solution allows a transparent scaling of performances on different devices. moreover, if the machine is equipped with more than one gpu, it is also possible to distribute the workload by launching the kernel on each gpu. right side. memory hierarchy: in cuda there are many different memories with different scopes. each thread has two different kind of private memory: registers and local memories. threads belonging to the same block can communicate through the shared memory, which has low access latency. the global memory suffers from high access latencies but it is accessible to all threads and it is cached since the introduction of the fermi architecture. also the texture and the constant memory are equipped with a cache as well, and all threads can read from these two memories given the peculiar features of cuda architecture, gpu programmers should be able to optimize both threads partitioning and memory usage, as well as to redesign the algorithm with appropriate kernels, in order to fully leverage the computational power of these multi-core devices. for instance, in the implementation of lassie, the shared memory is not used because several blocks are exploited to solve the system of odes, and threads do not communicate data with each other. moreover, the data structures employed by lassie are larger than the total size of the shared memory, thus preventing the possibility to exploit it. in what follows, we show how gpu programming and cuda features—including built-in support for vector types, which extend the standard c data types to vector—have been exploited to optimize the execution workflow of lassie. the paper is structured as follows. in the next section we briefly introduce the formalism of reaction-based models and provide a general description of lassie’s implementation. then, we discuss the computational performance of lassie, showing the speed-up it achieves with respect to lsoda for the simulation of reaction-based models of different sizes. we also analyze how the number of reactions and the number of species affect the performances of lassie. we conclude the work with some final remarks about cuda’s architecture and lassie, proposing future improvements of the simulator. lassie is available on the github repository https://github.com/aresio/lassie. methods lassie is designed to be a “black-box” deterministic simulator, created to be easily used without any gpu programming or odes modeling skills. in this section we describe how lassie allows to perform deterministic simulations of large-scale biochemical models, distributing all required calculations on the cores of the gpu. it is worth noting that the parallelization strategy exploited by lassie represents one of the novelties of this work, and allowed to achieve the remarkable performance results presented in the next sections. in particular, lassie has been developed to solve systems of coupled odes specified in the form dxdt=f, where x≡x represents the vector of concentration values at time t of all chemical species occurring in the system. reaction-based models and odes generation reaction-based modeling is a mechanistic, quantitative and parametric formalism to describe and simulate networks of biochemical reactions, which was exploited to analyze different signal transduction pathways. a reaction-based model is defined by specifying the set of n molecular species {s …,s n} and the set of m biochemical reactions {r …,r m} which appear in the cellular process under investigation. a generic reaction is described as follows: ri:∑j=1naijsj→ki∑j=1nbijsj,i= …,m, where a ij, bij∈ℕ are the stoichiometric coefficients and ki∈ℝ+ is the kinetic constant associated with r i. the set of reactions {r …,r m} can be written compactly in the matrix-vector form as→kbs, where s=[s 1⋯s n]t is the n-dimensional column vector of molecular species, k=[k 1⋯k m]t is the m-dimensional column vector of kinetic constants, and a,b∈n m×n are the so-called stoichiometric matrices whose elements i,j and i,j correspond to the stoichiometric coefficients a ij and b ij of the reactants and the products of all reactions, respectively. since a reaction simultaneously involving more than two reactants has a probability to take place almost equal to zero, here we consider only first and second-order reactions. for this reason, the matrices a and b are sparse. given an arbitrary reaction-based model and assuming the law of mass-action, it is possible to derive the corresponding system of coupled odes that describes the variation in time of the species concentrations. specifically, by denoting the concentration of species s j at time t as x j, where x j∈r ≥ for j= …,n, the system of coupled odes can be obtained as follows: dxdt=t, where x is the n-dimensional vector of concentration values at time t, the symbol ∘ denotes the entry-by-entry matrix multiplication, and x a denotes the vector-matrix exponentiation form. formally, x a is a m-dimensional vector whose i-th component is given by x1ai1⋯xnain, for i= …,m. we highlight that each ode appearing in eq. is a polynomial function, consisting in at least one monomial that is associated with a specific kinetic constant. data structures and cuda memory usage given a reaction-based model as input, lassie automatically generates the systems of odes according to eq. and encodes the matrices a and h=t as two arrays of short cuda vector types, named v a and v h, respectively. cuda vector types are multi-dimensional data ranging from to components, addressed by.x,.y,.z, and.w. since the matrices a and h are sparse, lassie uses compressed data structures created by removing all zero elements from a and h, in order to save memory and avoid unnecessary readings from the global memory. namely, let h ji be the element of h at row j and column i, and a ij the element of a at row i and column j, for i= …,m and j= …,n. for each non-zero element of h, we store into the.x and.y components of v h the values j and i, respectively; the.z component of v h is used to store the element h ji, while the.w component stores the index of the kinetic constant associated with that monomial. similarly, for each non-zero element of a, the.x and.y components of v a contain the values i and j, respectively. the value a ij is stored into the.z component of v a, while the.w component is left unused. note that we exploited the short cuda vector type rather than the short cuda vector type, because the former is 8-aligned and requires a single instruction to fetch a whole entry, while the latter is 2-aligned and thus takes three memory operations to read each entry. in order to parse these arrays inside the gpu, we use two additional arrays of short cuda vector types, named o h and o a, which store the offsets used to correctly read the entries of the v h and v a structures, respectively. the.x and.y components of each row of o h contain, respectively, the first index and the last index to access the v h structure. each thread uses its own pair of indexes to read the rows of the v h structure between the first index and the last one. similarly, o a stores the indexes that allow to correctly access the v a structure. finally, the values of the kinetic constants are stored into an array of type double, named k. figure shows an example of the matrix encoding used in lassie. fig. example of matrix encoding to automatically generate an ode using lassie. all terms of the polynomial function describing the ode of species x given at the top of the figure are encoded in the components of the data structures o h, v h, o a, v a and k, as detailed hereby. notice that only the data structures components with solid borders are used to automatically generate the ode; the various terms appearing in the ode are represented with corresponding colors in the data structure components. matrix encoding starts from matrix o h. each thread j, for j= …,n− reads the values stored in the.x and.y components of o h. in this example, we consider species x that corresponds to thread each thread fetches the values in v h, starting from the row indicated by the value stored in the.x component of o h, up to the row corresponding to the value stored in the.y component. in this example, thread in matrix o h reads the values contained in the first two rows—i.e., rows and 1—in matrix v h. each row of v h encodes a monomial of an ode: the.x component is not used; the.y components indicate the row numbers of the o a structure that each thread must read; the.z components indicate the sign and the coefficient of the monomial; the.w components indicate the positions of the array k containing the values of the kinetic constants corresponding to the reactions that the threads are parsing. in this example, the.z and.w components of v h allows to derive the coefficients −1k and +1k for the first and the second term of the ode, respectively. afterwards, as in the case of o h, each thread fetches the values in v a, starting from the row indicated by the value stored in the.x component of o a, up to the row corresponding to the value stored in the.y component of o a. the values stored in the.y and.z components of v a correspond to the indexes of the species and the stoichiometric coefficients, respectively, while the.x and.w components of v a are left unused. in this example, row in matrix o a reads the values stored in rows and of matrix v a, generating the factors in the first term of the ode, while row in matrix o a reads the values stored in row of matrix v a, generating the factor in the second term of the ode. therefore, in this example, the matrix encoding overall generates the ode of species x consisting in the sum of two polynomial terms: −k 111+k 21 thanks to these cuda structures, we obtain a twofold performance improvement: at the instruction level, a single instruction is enough to either load or store a multi-word vector. so doing, the total instruction latency for a particular memory transaction is lower and also the bytes per instruction ratio is higher; at the memory controller level, by using vector types a transfer request from a warp has a larger net memory throughput per transaction, yielding a higher bytes per transaction ratio. with a fewer number of transfer requests, the memory controller is able to reduce contentions producing a higher overall memory bandwidth utilization. the only limitation due to short data type is that indices are limited to 22×8− which means that lassie cannot simulate systems larger than chemical species and reactions. execution workflow and cuda kernels once that the system of odes is generated by reading the input files and appropriately stored according to the cuda vector types, lassie solves it by automatically switching between the runge-kutta-fehlberg method in the absence of stiffness, and the backward differentiation formulae methods in presence of stiffness. the integration of the systems of odes is carried out from an initial time instant t up to a given maximum simulation time t max. in order to reproduce the dynamics of the cellular process described by the odes, the concentration values of the molecular species appearing in the reaction-based model are saved at specified time steps within the interval [t t max]. lassie’s workflow consists in distinct phases, as represented in fig. note that phases p p and p are executed by the host, while p p and p are executed by the device. overall, phases p p and p rely on different lightweight kernels, which were specifically developed to fully leverage the parallel architecture of the gpu for the implementation of the aforementioned numerical integration methods. we describe hereafter the main design and implementation choices of each phase and their related cuda kernels, which result in a novel parallelization strategy with respect to state-of-the-art methodologies. fig. simplified scheme of lassie workflow. the data structures used to encode the system of odes are generated in phase p in phase p if the current simulation time t corresponds to a specified sampling time instant, then the current concentration values of all molecular species are saved; otherwise, the execution proceeds to the next phase. in phase p each thread derives and solves the corresponding ode by exploiting the rkf method, while in phase p the rkf solutions are verified: if the rkf solutions are rejected, then the integration step-size dt is reduced and phase p is executed again; if rkf solutions are rejected but the integration step-size dt is too small, then phase p is executed and the system of odes is solved using the bdf methods; if the rkf solutions are accepted, the termination criterion is verified during phase p phase p 1. it implements the generation of all data structures used to encode the odes, as described in the previous section. this phase is executed on the host. phase p 2. it is used to sample and save the system dynamics and it is implemented by means of a single cuda kernel. in particular, if the current simulation time t corresponds to one of the specified sampling time instants, lassie saves the concentration values of all molecular species into an array defined on the gpu. otherwise, the execution proceeds to the next phase. phase p 3. it implements the rkf method, an explicit integration algorithm with variable step-size used by each thread j to solve the j-th ode, for j= …,n− this phase is implemented as cuda kernels. during this phase, two different approximated states u and w of the state x of the system are generated at each step, thanks to the evaluation of six supplementary values l …,l. to evaluate the accuracy of u and w at the current step-size dt, lassie exploits a user-defined vector tolerance ε∈r n, and two additional arrays, e r,δ∈r n, defined as follows: er=|w−u|dt,δ= εer if e r j≤ε j for all j= …,n, then u is accepted as new state of the system, that is, x=u; otherwise, the solutions u and w are rejected and recalculated by using a new step-size. the new step-size is computed as d t=d t· min{δ …,δ n}, being δ …,δ n the components of vector δ. overall, phase p is implemented by means of the following kernels: kernel k 2: used to evaluate each ode at the current state x of the system; kernels k – k 8: each thread j, for j= …,n− computes the components l 1j,…,l 6j of l …,l by invoking kernel k 2; kernel k 9: each thread j, for j= …,n− computes the components w j and u j of the approximated states u and w, respectively; kernel k 10: each thread j, for j= …,n− calculates the components e r j and δ j of e r and δ, respectively. phase p 4. it is used to verify the rkf solutions calculated during phase p and, accordingly, to choose the next phase to be executed: if the solutions are rejected and the new step-size dt is acceptable, phase p is executed again exploiting a smaller step-size dt; if the solutions are rejected and the new step-size dt becomes too small, lassie executes phase p 5; if all solutions do not violate the specified rkf-tolerance vector ε, then lassie executes phase p note that point implicitly states that the system of odes is considered to be stiff, so that lassie automatically switches to phase p where bdf methods are used for the numerical integration. phase p is executed on the host. phase p 5. it implements the bdf methods, the most widely used implicit multi-step numerical integration algorithms. lassie switches to this phase if and only if the rkf solutions u and w evaluated during phase p are rejected, and the rkf step-size dt becomes smaller than ε s. the general formula for a bdf can be written as ∑i=0qαix=dtβ0f), where the coefficients α i and β are chosen according to the order q of bdf, and dt is user-defined. note that, for q> the absolute stability region of the resulting bdf methods is too small and such bdfs are numerical unstable. therefore, bdfs with an order q greater than are not used. since each bdf is an implicit method, at each time step it requires the solution of a nonlinear system of equations, which can be solved by using the iterative newton—raphson method. this algorithm allows to find successively better approximations z of the zeros of a real-valued function f= by using the derivative of f, and it is repeated until a sufficiently accurate value is reached. this idea can be extended to a system of nonlinear equations, by using the jacobian matrix j) of f), which is the matrix of all first-order partial derivatives. since the evaluation of the jacobian matrix at each iteration is computationally expensive, lassie actually exploits: a modified newton—raphson method; the lu factorization method. during phase p the newton-raphson method is iterated until a user-defined maximum number of iterations m a x it is reached, or a sufficiently accurate value is achieved. overall, phase p is implemented by means of the following kernels: kernel k 11: each thread j, for j= …,n− derives the j-th row of the jacobian matrix and evaluates it on the current state of the system x; kernel k 12: the jacobian matrix is transposed in order to exploit the lu factorization method; kernels k – k 18: based on the order q of the bdf, lassie invokes one of these kernels to calculate the known terms of the linear system; kernels k – k 24: each kernel k, q= …, performs the calculations of the q-th order bdf; kernel k 25: it updates the iteration vector needed to execute the newton-raphson method. phase p 6. it is used to verify the termination criterion: if the maximum time t max is reached, then the simulation ends. on the contrary, the execution iterates from phase p this phase is executed on the host. all the temporary results computed by lassie are stored on the gpu, since data transfers between the host and the device are very time consuming. for the same reason, the output data are transferred to the host as soon as the whole simulation is completed.  in this section we compare the computational performance of lassie against lsoda, which is generally considered one of the best numerical integration algorithms for deterministic simulations of biological systems, thanks to its capability of dealing with stiff and non-stiff systems. in particular, we exploited the lsoda implementation provided by scipy library, written in c language. lassie was run on a machine with a gpu nvidia geforce titan gtx, based on the kepler architecture and equipped with 2× streaming multiprocessors for a total of cores and a theoretical peak processing power of tflops in double precision. instead, lsoda was run on galileo, a supercomputer created by the italian consortium cineca. galileo consists of compute nodes, each one equipped with cpus octa-core intel xeon haswell e5- v for a total of cores, and gb of ram. each cpu is capable of about gflops in double precision. in our tests, we exploited one node with gb of ram distributed over cores. the computational performance was evaluated by simulating a set of synthetic reaction-based models of increasing size, that is, having a number of reactions and species m×n arbitrarily chosen in the range from 64× to 8192× the models were generated considering the methodology used in, which was modified in order to randomly sample the initial concentration of each species with a uniform distribution in the range [ ), and the kinetic constant of each reaction with a logarithmic distribution in the range [10− ). for each model size m×n, we generated and simulated different synthetic reaction-based models to the aim of measuring the average running time of both lassie and lsoda. the simulation of each reaction-based model was performed multiple times, using different settings for the sampling of the time-series. specifically, in each repetition, we saved either, or samples of the system dynamics of all chemical species, at regular intervals. all simulations were halted at time t max=. all simulations were executed—independently from the size of the model and the number of samples saved—by setting the following parameters of lassie: tolerance of rkf method ε j = 10− j= …,n; first–order bdf method; bdf integration step d t=; tolerance of newton-raphson method ε nr = 10−6; maximum number of iterations allowed during each call of the newton-raphson method m a x it=104; initial integration step of rkf method equal to 10−3; tolerance value to switch between rkf and backward euler methods ε s=10− the following parameters of lsoda were used to run the simulations: relative tolerance equal to 10−6; absolute tolerance equal to 10−12; maximum number of internal steps equal to fig. speed-up values achieved by lassie with respect to lsoda for the simulation of synthetic models of increasing size, having a number of reactions and of species m×n and characterized by an increasing number of sampling time instants of the system dynamics. when the value of the speed-up is greater than one, lassie is faster than lsoda and vice versa m×n lsoda lassie speed-up lsoda lassie speed-up lsoda lassie speed-up lsoda lassie speed-up lsoda lassie speed-up na na na na na na na na na na *maximum speed-up value fig. comparison between the average running time required by lassie and lsoda to simulate instances of models characterized by reactions and species, reactions and species, reactions and species, saving different numbers of sampling time instants of the dynamics. note that the y-axes are in logarithmic scale as an additional test, we investigated whether the relationship between the number of reactions and the number of species could affect the overall performances of lassie. as the number of chemical species corresponds to the number of odes, the length of each ode is roughly proportional to the number of reactions. since gpus have a lower clock frequency than cpus, each gpu core is slower than the cpu core to perform a single instruction for this reason, in order to obtain the highest performances, the calculations on the gpu should be spread across threads as much as possible, while the number of operations performed by each thread should be reduced. indeed, as reported in table and shown in fig. when the number of chemical species involved in a model is greater than the number of reactions, lassie achieves better performances than those obtained in the case of models with a number of chemical species smaller than the number of reactions. for instance, considering the models with m×n equal to 171× the running time of lassie is smaller than in the case of the models with size 512× irrespective of the number of samples of the system dynamics, thanks to the higher number of threads that are concurrently launched on the gpu in the first case. this is in general valid in all cases with the exception of the models characterized by chemical species with and samples of the system dynamics. here, the average running time of lassie is greater than in the case of models with reactions, since the required number of accesses to the high-latency global memory of the gpu impairs the performances of the simulations. fig. running time of lassie for the simulation of synthetic models of increasing size, having a number of reactions and of species m×n, with m≠n, and characterized by an increasing number of sampling time instants of the system dynamics m×n lassie lassie lassie lassie lassie in order to assess the scalability of lassie, and of cuda applications in general, we executed additional tests on different gpus. figure shows a comparison of lassie’s performance using three different gpu models: a notebook video card, the nvidia geforce titan z used throughout the paper, and a tesla-class gpu. to compare the speed-up provided by these gpus we generated different synthetic models and calculated the average running time. fig. comparison of the average running times for the simulation of synthetic models characterized by three different sizes, executed with different gpus: a notebook gpu nvidia geforce 960m; a nvidia geforce titan z; a tesla-class gpu nvidia k20c our results highlight the importance of two distinct factors on lassie’s performances: the gpu’s clock frequency and the amount of available resources. as a matter of fact, despite the lower amount of cuda cores, the geforce 960m turns out to be competitive on models of moderately large size thanks to its higher clock rate, with respect to the titan z and the k20c. when the odes largely outnumber the available cores, the geforce 960m is no longer competitive. this is an example of transparent scalability of cuda applications: the threads are automatically distributed over the available cores, improving the overall performances, without any user intervention. moreover, as described in the background section, threads are organized in blocks that are scheduled on the available multi-processors. thanks to this characteristic, when the overall number of threads outnumbers the available cores, cuda automatically creates a queue of blocks that are scheduled on the streaming multi-processors as soon as they become available for computation. thus, lassie can, in principle, simulate any model on any gpu, as long as there is enough memory to store the data structures. the tesla k20c is characterized by a large amount of cores that, in the case of 1024× models, are fully exploited only during the simulation of the stiff parts of the dynamics. for the remaining parts of the simulation, half of its cores are actually used for computation with a slower clock rate with respect to the clock rate of the geforce gpus. moreover, tesla cards exploit error correcting codes on memories, ensuring additional checks of correctness to the data against potential corruption from electrical or magnetic interference, at the price of a significant overhead. the ecc was enabled during all tests, partly explaining the reduced performance of the tesla k20c on very large-scale models with respect to the titan z. we assessed the accuracy of lassie by simulating the dynamics of the model of the ras/camp/pka signaling pathway in yeast presented in, and comparing the outcome of lassie with the result of the simulation performed with lsoda. we also investigated the influence of lassie parameters on the running times and quality of the simulated solutions, by exploiting a model representing a chain of isomerizations. the accuracy results—which show an identical dynamics with respect to lsoda using default settings—are presented in the additional file as a final remark, we highlight that a fair comparison of gpus and cpus is a difficult task, in general, due to their deep architectural differences. the theoretical peak performances of both architectures are difficult to achieve: indeed, developers must implement code to the aim of maximizing the parallelism and the occupancy of the multi-processors, adhering as much as possible to the underlying simd computational model in the case of the gpu and exploiting vector instructions in the case of the cpu. however, gpus allow the temporary divergence of the execution flow of threads, that is, a part of the threads can execute different portions of the code. when this situation occurs, some threads get stalled waiting for reconvergence. this mechanism provides the programmer with a certain degree of freedom to abandon the simd paradigm, but at the same time it can potentially lead to the complete serialization of the execution affecting the overall performances. hence, conditional branches should be avoided as much as possible. we also highlight that the usage of registers and shared memory influences the occupancy of the gpu, as these resources are scarce on each streaming multiprocessor. all these circumstances can prevent the achievement of the peak computational power of a gpu. to this aim, we developed kernels that maximize the parallelism and the occupancy of the multi-processors avoiding threads divergence as much as possible. moreover, we optimized data structures to store the matrices a and h that encode the system of odes, and cuda vector types that allow to increase the memory throughput and to reduce the number of memory accesses, all precautions that explain the performance boost achieved with lassie.  in this work we presented lassie, a gpu-powered simulator of large-scale biochemical systems based on mass-action kinetics. lassie is a “black-box” simulator able to automatically convert reaction-based models of biological systems into the corresponding systems of odes. reaction-based models defined according to the law of mass-action do not hinge upon the use of any approximate kinetics functions, which are frequently used in systems biology for the definition of mathematical models based on differential equations. although michaelis-menten kinetics or hill functions can be useful in biological modeling, they rely on chemical assumptions that are valid only in certain conditions. therefore, the reason why we rely on mass-action based models is manifold. on the one hand, since the biological function and biochemical kinetics of all molecular species and all reactions appearing in the model are not approximated nor lumped together in any way, they can be analyzed independently from each other. as a consequence, this allows to determine the influence of every single species and reaction on the overall functioning of the system. on the other hand, the law of mass-action allows to derive a first order ode for each species appearing in the model: it is worth noting that such ode is a polynomial function that describes how the concentration of that species changes in time, according to all the reactions where it appears either as reactant or product. the presence of polynomial functions simplifies the symbolic derivation that is needed to calculate the jacobian matrix associated with the odes and exploited by the bdf. in addition, as described in the methods section, polynomials can be efficiently encoded in the memory and parsed gpu-side. as a result, all gpu threads can perform the same task, strongly reducing warps’ divergence and the consequent stalling of threads due to serialization, a circumstance that would instead happen if each thread calculated an ode characterized by an arbitrary kinetics. in order to solve systems of odes characterized by stiffness, lassie automatically switches between the rkf and the bdf integration methods. lassie’s execution flow is partitioned into cuda kernels, overall distributing the calculations over the available cores in order to fully exploit the massive parallel capabilities of modern gpus, therefore achieving a relevant reduction of the running time in case of large-scale models. in order to assess the computational performance of lassie, we performed a set of simulation tests using synthetic reaction-based models of increasing size, and we compared lassie’s running time with respect to the lsoda numerical integration algorithm implemented in the scipy library. the break-even between the performances of lassie and lsoda was observed when both the numbers of reactions and chemical species is in between and this result indicates that, for biological systems consisting in more than reactions and species, the gpu-powered simulator becomes more convenient than the lsoda algorithm running on cpu. indeed, in the case of large-scale models, characterized by reactions and species, we obtained a considerable 92× speed-up. moreover, thanks to its smaller memory footprint with respect to lsoda, lassie allows the simulation of even larger models, taking just an average of s to simulate models characterized by reactions and species. on the contrary, lsoda did not allow the simulation of models of this size on the computer we used for the tests, as it crashed because of its very high memory footprint. we also highlight that copasi, one of the most used software in systems biology, requires in general longer execution times with respect to the scipy implementation of lsoda exploited in this work. in addition, copasi fails when trying to simulate large-scale models. we provide an example of such model—characterized by reactions and species—as sbml file in the github repository. bdfs are the most used integration algorithms to solve systems of odes in case of stiffness. the first–order bdf is a single-step implicit integration method, meaning that the next state of the system depends only on the current state of the system. higher–order bdfs are multi-step methods, so that the next state of the system relies on multiple previous states of the system. this implies that the integration step-size should be the same for all previous states to ensure the correctness of the solution. for this reason, lsoda uses the multi-step adams methods as explicit methods in addition to bdfs. conversely, lassie uses the rkf method, which is a single-step explicit algorithm with variable step-size, and the backward euler method. other single-step implicit methods belonging to the family of runge-kutta methods exist, the most known being the families of lobatto and radau methods. these methods have been proven to be suitable for stiff systems, thanks to their accuracy and stability. as a future development of lassie, we will investigate the feasibility and efficacy of replacing the backward euler and, more generally, the bdfs with implicit runge-kutta methods, in particular lobatto and radau methods. in order to fully exploit the cuda architecture, the memory hierarchy must be exploited as much as possible. because of the peculiar sequential structure of both explicit and implicit integration algorithms, lassie’s kernel are lightweight and rarely reuse any variables. for this reason, the current implementation only leverages the global memory and registers to manipulate the mutable data. the shared memory has not been exploited in any way, leaving room for potential future improvements of performances. however, on the gpus where the l cache and the shared memory share the same resources, cuda allows to express a preference to assign a larger amount of memory to the caching mechanisms. this functionality is enabled by default in lassie using the cuda cudafuncsetcacheconfig primitive, executed with the cudafunccachepreferl argument. also the constant memory has not been used, since all data structures are larger than the total size of this memory. in a future release, we plan to leverage these memories, for instance to store the array of the kinetic constants and the structures containing the offsets used to correctly decode the odes. lassie currently exploits only a single gpu, even on multi-gpu systems; as a future improvement of this work we plan to extend it in order to support multi-gpu systems, to further increase the size of the models to be simulated. finally, an additional goal in the development of lassie is to integrate and accelerate the investigation of rule-based models. in the next future, we plan to develop a set of tools that will leverage the functionalities offered by rule-based modeling frameworks, to convert a rule-based model into a set of reactions. although rule-based tools already provide internal simulation methods, lassie can represent a valuable alternative for large-scale models, enabling the investigation of more detailed biological systems, paving the way to potential new discoveries in systems biology. lassie will be also integrated in cosys, a free web-based platform for systems biology investigation available at http://www.sysbio.it/cosys. endnote the advances in gpu’s technology will progressively reduce this gap. in middle 2016—with the introduction of the novel pascal architecture and the nm finfet manufacturing process—nvidia presented a gpu with a clock frequency of ghz that, theoretically, is expected to double lassie’s performances. additional files additional file lassie graphical user interface. additional file lassie input files and command line arguments. additional file implementation of cuda kernels for lassie execution workflow. additional file simulation accuracy of lassie. abbreviations bdfbackward differentiation formulae cpucentral processing unit cudacompute unified device architecture eccerror correcting codes gpgpugeneral-purpose gpu gpugraphics processing unit lassielarge-scale simulator odeordinary differential equation rkfrunge-kutta-fehlberg simdsingle instruction multiple data ssastochastic simulation algorithm 
notably, thanks to its smaller memory footprint, lassie is able to perform fast simulations of even larger models, whereby the tested cpu-implementation of lsoda failed to reach termination. mathematical modeling and in silico analysis are widely acknowledged as complementary tools to biological laboratory methods, to achieve a thorough understanding of emergent behaviors of cellular processes in both physiological and perturbed conditions. given a reaction-based model of a cellular process, lassie automatically generates the corresponding system of ordinary differential equations, assuming mass-action kinetics. though, the simulation of large-scale models—consisting in hundreds or thousands of reactions and molecular species—can rapidly overtake the capabilities of central processing units. lassie is a “black-box” gpu-accelerated deterministic simulator, specifically designed for large-scale models and not requiring any expertise in mathematical modeling, simulation algorithms or gpu programming. lassie is therefore expected to make an important breakthrough in systems biology applications, for the execution of faster and in-depth computational analyses of large-scale models of complex biological systems. the computational performance of lassie are assessed using a set of randomly generated synthetic reaction-based models of increasing size, ranging from to reactions and species, and compared to a cpu-implementation of the lsoda numerical integration algorithm. by virtue of this implementation, lassie achieves up to 92× speed-up with respect to lsoda, therefore reducing the running time from approximately month down to h to simulate models consisting in, for instance, four thousands of reactions and species. lassie adopts a novel fine-grained parallelization strategy to distribute on the gpu cores all the calculations required to solve the system of odes.
 mathematical modeling and in silico analysis are widely acknowledged as complementary tools to biological laboratory methods, to achieve a thorough understanding of emergent behaviors of cellular processes in both physiological and perturbed conditions. though, the simulation of large-scale models—consisting in hundreds or thousands of reactions and molecular species—can rapidly overtake the capabilities of central processing units. the purpose of this work is to exploit alternative high-performance computing solutions, such as graphics processing units, to allow the investigation of these models at reduced computational costs.  lassie is a “black-box” gpu-accelerated deterministic simulator, specifically designed for large-scale models and not requiring any expertise in mathematical modeling, simulation algorithms or gpu programming. given a reaction-based model of a cellular process, lassie automatically generates the corresponding system of ordinary differential equations, assuming mass-action kinetics. the numerical solution of the odes is obtained by automatically switching between the runge-kutta-fehlberg method in the absence of stiffness, and the backward differentiation formulae of first order in presence of stiffness. the computational performance of lassie are assessed using a set of randomly generated synthetic reaction-based models of increasing size, ranging from to reactions and species, and compared to a cpu-implementation of the lsoda numerical integration algorithm.  lassie adopts a novel fine-grained parallelization strategy to distribute on the gpu cores all the calculations required to solve the system of odes. by virtue of this implementation, lassie achieves up to 92× speed-up with respect to lsoda, therefore reducing the running time from approximately month down to h to simulate models consisting in, for instance, four thousands of reactions and species. notably, thanks to its smaller memory footprint, lassie is able to perform fast simulations of even larger models, whereby the tested cpu-implementation of lsoda failed to reach termination. lassie is therefore expected to make an important breakthrough in systems biology applications, for the execution of faster and in-depth computational analyses of large-scale models of complex biological systems. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords graphics processing unitgpu computingreaction-based modeldeterministic simulationnumerical integration methodlsodanvidia cudafine-grained parallelizationsystems biologyrule-based modelissue-copyright-statement© the author 2017 