Mex code generator

From BioNetWiki

writeMexfile() is a new BioNetGen action that generates a model-specific ODE integrator in the C language. The integration is implemented using the Sundials CVode library, while data input and output is managed through a MATLAB Mex interface that wraps around the integrator. The CVode-MEX implementation is both efficient and flexible. After compilation, the model integrator is useful for tasks that require many repeated simulations of the model with various parameters, times, and initial conditions (such as parameter estimation, model sensitivity, model checking, etc). The user must have access to a recent version of MATLAB in order to compile and run the model integrator! Additionally, the freely available Sundials CVODE library (v2.6.0) must be installed.

writeMexfile is a BETA feature. If you encounter problems, please send us a detailed description (better yet, send the console output, BNGL model, and output files). Your feedback is welcome.

Instructions

Before creating CVode-Mex code for your BNG model, be sure to generate, debug and validate the model reaction network. Then call the "writeMexfile" action command:

   writeMexfile({});

There are a number of arguments that writeMexfile accepts, including:

argument	description	default value
prefix	file prefix	[model_name]
atol	absolute scalar tolerance in CVODE	1e-6
rtol	relative scalar tolerance in CVODE	1e-8
t_start	default start time of the simulation	0.0
t_end	default end time of the simulation	10.0
n_steps	number of observable output steps, not including the initial timepoint	10
stiff	select the stiff CVODE solver	true=1
sparse	select the sparse jacobian approximation in CVODE (GMRES)	false=0
max_step	max step size in CVODE	0.0, i.e. no maximum
max_num_steps	max number of CVODE integration steps	500
max_err_test_fails	max number of CVODE error test failures	7
max_conv_fails	max number of CVODE convergence failures	10

Note that t_start, t_end, and n_steps can be changed at runtime. These options just set the defaults.

For example:

   writeMexfile({atol=>1e-4,t_start=>0,t_end=>1000,n_steps=>100,sparse=>1})
   

Will generate a CVode-Mex model integator with absolute tolerance of 1e-4 that integrates (by default) from t=0 to t=1000 using the GMRES sparse method.

writeMexfile creates two files in the current directory:

  [prefix]_mex.c
  [prefix].m
  

The first file contains source code that defines the model integrator and Mex-wrapper. The source must be compiled using the "mex" function in MATLAB. The second file is a matlab script that sets default parameter values and initial conditions, then calls the CVode model, and finally generates a plot of observable trajectories.

To compile the mex files (GNU/Linux): You must have the Sundials CVode library installed in a place where MATLAB can find it. The model integrator was tested using CVode v2.6.0 and I cannot guarantee it will work with other versions (including the version packaged with BNG!). Open Matlab and navigate to the directory with the mex file. Call the "mex" compiler as follows:

  mex -lsundials_nvecserial -lsundials_cvode -lm [model_name]_mex.c

Matlab may issue a version warning, but the file should compile without any fatal errors. If Matlab can't find the Sundials libraries or header files, tell Matlab where to look using the -L and -I switches.

On the Mac OS/X and Windows platforms, follow John Sekar's instructions for Compiling Mex files (pdf).

After successful compilation, the Mex file can be called directly from the MATLAB console using the following syntax:

  [err_flag, species, obsrv] = [prefix]_mex( params, init_species, timepoints );

The input arguments are a column vector of model parameters, a row vector of initial species counts, and a row vector of timepoints. The output arguments are an error flag (0 if the model integrated successfully), an array of species counts (rows correspond to timepoints), and an array of observable counts.

For heavy computational tasks, it's best to call the model integrator directly. But defining the parameters an initial conditions can be unwieldy at times. The companion script [prefix].m simplifies this process! The script contains the default parameter values and mathematical expressions defining initial species counts in terms of parameters. This eliminates the tedious task of copying parameters and mathematical expressions from the BNGL model into MATLAB. To simulate the model using default parameters and initial conditions, simply call the script with empty arguments:

  [err_flag, species, obsrv] = [prefix]( [], [], [] );
  

The script will call the model integrator and create a figure containing observable trajectories. The parameters and timepoints may be changed by passing a column vector of parameters and a row vector of timepoints:

  [err_flag, species, obsrv] = [prefix](params, [], timepoints );
  

The initial concentrations are calculated from the new parameters using the expressions defined in the BNGL model. If different initial conditions are desired, pass a row vector with the new initial counts:

  [err_flag, species, obsrv] = [prefix]( params, init_sp, timepoints );
  

The script outputs the species and observble trajectories, permitting further analysis and custom plots.

Multi-stage simulations are simple: Just call the script again using the species vector obtained from the last timepoint of the previous stage:

  [err_flag1, sp1, obsrv1] = [prefix]( params1, init_sp1, [0:1:100] );
  [err_flag2, sp2, obsrv2] = [prefix]( params2, sp1(end,:), [100:1:200] );

Notes

1. The BNGL parameter block includes free parameters and derived parameters. A derived parameter depends on other parameters and is essentially "fixed" once the free parameters are set. writeMexfile() differentiates between free and derived parameters. The CVode-Mex model integrator only accepts free parameters as input. Derived parameters are calculated at run-time using the same mathematical expressions defined in the BNGL model file. If the user wishes to change a derived parameter in the model, then the mathematical expression in the BNGL file should be replaced with a numerical value.
2. The CVode-Mex code generator supports elementary ratelaws, global functions, and local functions. Support for the "special" ratelaws (Sat, MM, Hill) is not yet available (but these ratelaws can be implemented using global functions).
3. Ratelaws are always interpretted as microscopic rates (i.e. per reactant set). This conforms to the BioNetGen convention for elementary ratelaws, but contrasts the NFsim convention for global functions. Be careful defining non-elementary ratelaws. The common Michaelis-Mentin ratelaw, in microscopic terms, should be written as:

    MM() = kcat*Etot/(km + S_total)

   where S_total is an observable that matches the total quantity of substrate. S_total does NOT appear in the numerator because the function defines the "per molecule of substrate" rate. Thus, the overall rate is equal to the product of this function and the substrate cooncentration (and the S_total magically reappears in the numerator, yieling the conventional MM ratelaw).
4. Future implementations will include an option to generate MATLAB independent code.