function [t,n] = gma(tspan,pars)
NAV = pars{1};
S0 = pars{8};
options = odeset('refine',16);
[t,S] = ode15s(@compdS,tspan,S0,options,pars); % solution of gma
n = NAV*S*1e-9; % number of molecules predicted by gma
% compute derivative for gma
function dS = compdS(t,S,pars)
k = pars{2};
M = pars{4};
p = pars{5};
l = pars{6};
nu = pars{9};
for u=1:M    
    P(u) = prod(S(p{u}(:)).^l{u}(:));
end
dS = ( k.*P*nu )';