Tutorial¶
For this tutorial, we will fit the parameters to the following ODE system:
Note
The Python script and datafile for this tutorial is available at simpleKinetics_example.py and simpleKinetics_data.mat, respectively.
The data for simpleKinetics_data.mat were generated using the parameters:
Suppose we have no idea what the k parameters could be. We establish a very generous uniform prior on the parameters:
- k1 bounded at [0,5]
- k2 bounded at [0,7]
Initialization of an APTModel Object¶
To begin, we need to import the APT-MCMC python package, set the locations of our installed Sundials and TRNG libraries, as well as initialize the APTmodel object. This can be done by:
>> import aptmodel as apt
>> model = apt.aptmodel() # Initialize model
>> model.set_include_path('/usr/local/include') # Set include dir of CVODE and TRNG libraries
>> model.set_lib_path('/usr/local/lib') # Set lib dir of CVODE and TRNG libraries
Creating parameters¶
Next, we need to tell our APTmodel object what the parameters are. The structure to do this is model.add_par(‘parname’, lower_bound, upper_bound).
Note
Currently, only uniform priors are supported. However, TRNG supports a variety of probability distributions. Consult the TRNG documentation to see an exhaustive list of PDFs supported and the APTModel documentation for how to change the prior distribution type in the C++ code manually.
For our system, we do the following:
>> model.add_par('k1',0,5) # define a parameter, k1, whose searchspace lies between (0,5)
>> model.add_par('k2',0,7) # define a parameter, k2, whose searchspace lies between (0,7)
>> model.showpar() # See the parameters and verify them before moving on
===aptmodel object, showing Parameters===
The following parameters and bounds are loaded:
Par 1: k1 , LB= 0, UB= 5
Par 2: k2 , LB= 0, UB= 7
>> model.output('simpleKinetics_example') # Set the name of the subfolder where the generated C++ code will live.
Note
The default output directory is ./Output. APTmodel WILL overwrite anything in the specified output directory. Therefore, it is not recommended to leave model.output at its default values. Each simulation should have its own directory.
Defining ODES¶
Before we define the ODEs, we have to give APTModel a list of the state names:
>> model.add_states('Ca','Cb') # 'Cc', 'Cd', 'Ce'.... etc
>> model.showstates() # See the states and verify them before moving on
===aptmodel object, showing States===
The following states are set:
State 1: Ca
State 2: Cb
Now that APTModel knows what the parameters are and the states are, we can enter in the ODE. The structure to do this is model.ODE[‘state_i’].set(‘C-based math expression’).
APTModel does not have the functionality to parse a Python math expression and convert it into C/C++-based math expression. Therefore, the responsibility is up to the user to write their ODEs in a suitable C/C++-based format. This is less difficult than it seems because the only significant difference is that Python uses x**y to raise x to the power of y, while C/C++ uses pow(x,y) .
Note
The following custom functions are built-in for user convenience:
- hill(a,b) = \(\frac{a}{a+b}\)
- hilln(a,n,b) = \(\frac{a^n}{b^n + a^n}\)
- heavy(z) = \(\frac{1}{10^{-10z}}\)
To define the equations in this example:
>> model.ODE['Ca'].set('5 - k1 * Ca')
>> model.ODE['Cb'].set('k1*Ca-k2*Cb')
Defining Experiments¶
Our APTmodel object, model, has now been fully defined. The next step is to add data for the simulation to fit. This is done by creating experiments, loading data to each experiment, and then finally adding the experiments into APTmodel.
For our tutorial, we have two experimental datasets that we want to fit to obtain the parameters for our model. We create the experiments, load data into it, and add them to our APTmodel object:
>> experiment1 = apt.experiment('Experiment1')
>> from scipy.io import loadmat
>> dataA = loadmat('simpleKinetics_data.mat')['dataA'] # dataA is nStates x nTime
>> # Create a APTdata object to store dataA in
>> aptdataA = apt.data(dataA,states=['Ca','Cb'], time=[0,2,4,6,8,12,24], isLog=False)
>> experiment1.add_data(aptdataA) # This function deepcopies aptdataA into experiment1 in case the user to redefines aptdata later.
>> model.add_experiment(experiment1)
Similarly for experiment 2, we:
>> experiment2 = apt.experiment('Experiment2')
>> dataB = loadmat('simpleKinetics_data.mat')['dataB'] # dataB is nStates x nTime
>> aptdataB = apt.data(dataB,states=['Ca','Cb'],time=[0, 2, 4, 6, 8, 12, 24],isLog=False)
>> experiment2.add_data(aptdataB)
>> experiment2.show(showall=True) # Check if experiment 2 was added correctly
===apt.experiment object: Experiment2===
nState x nTime aptdata matrix (2 x 7)
Preview (in linear space):
Ca Cb
t=0 0.00e+00 5.00e+00
t=2 8.41e+00 1.74e+00
t=4 1.43e+01 2.32e+00
t=6 1.84e+01 3.18e+00
t=8 2.13e+01 3.86e+00
t=12 2.47e+01 4.70e+00
t=24 2.76e+01 5.40e+00
==========================================================
>> model.add_experiment(experiment2)
In the above example, we could have also verified the data for experiment1 using experiment1.show(). The showall=True toggle overrides showing only the first and last few timepoints worth of data.
Note
Data for each state does not need to be loaded. APT-MCMC will only compare against available data. However, initial conditions for each state MUST be provided. If you are trying to fit against initial conditions, consult APTExperiment documentation.
Setting MCMC Options¶
Now we need to set APT-MCMC options by defining a separate options object:
>> options = apt.aptoptions(maxThr=8) # Create an options variable for model
>> options.filename = "simpleKinetics" # Set the name of the compiled executible
>> options.show() # See the values of all options. They can be set directly via options.OPTION_NAME = [value]
===aptoptions object===
Core Settings
|nTemp = 4, number of parallel ensembles to run
|nStep = 25, number of MCMC steps to take before attemping location swaps between temperatures
|nEnsemble = 100, number of samplers exploring each temperature simultaneously
|nSwap = 1e+04, number of swaps to attempt before ending simulation
|nInit = 1e+03, number of initial MC samples to take before starting simulation
.... (output truncated)
>> options.nSwap= 1e2 # Control simulation length by setting number of swaps to 100
Each option is briefly explained in the output. Consult APToptions documentation for more information.
Generating Simulation Files¶
All that is left is:
>> model.generate(options)
We now navigate to the output directory we set earlier in ./simpleKinetics_example/. From there, we run the following shell commands:
$ make
$ ./simpleKinetics
If all went well, compilation should have been successfull and you should see APT-MCMC running!
Loading Simulation Results¶
Once the simulation is done, we can load up the results in MATLAB or Python.
- MATLAB: MCMC_convert.m provides the function:
results_struct = MCMC_convert('filename.mcmc') - Python: APTMODEL package provides the function:
results_object = apt.openresults('filename.mcmc')