Brief Tutorial on Using DOTcvpSB (a Matlab Toolbox for by mll78346

VIEWS: 23 PAGES: 52

									                           B RIEF T UTORIAL ON U SING DOT CVP SB
          ( A M ATLAB T OOLBOX FOR DYNAMIC O PTIMIZATION IN S YSTEMS B IOLOGY )

                                                      Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga

                                                                        Process Engineering Group
                                                             Instituto de Investigaciones Marinas [IIM-CSIC]
                                                                 Spanish Council for Scientific Research
                                                                C/Eduardo Cabello 6, 36208 Vigo, Spain
                                                              e-mails: {thirmajer, ebalsa, julio}@iim.csic.es

                                                             All information about the toolbox is available at
                                                                     http://www.iim.csic.es/∼dotcvpsb/

                                                                   or partial information is available at
                                                                   http://reven.comli.com/index4.html




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga          ()            Brief Tutorial on Using DOTcvpSB      April 5, 2009   1 / 43
                                                                  Outline



   O UTLINE


          1    DOTcvpSB: Key Features
          2    Optimal Control Problems
          3    NLP and MINLP Solvers, Solution Scheme

          4    Single Optimization
          5    Multistart
          6    Sucessive Re-optimization
          7    Hybrid Strategy
          8    SBML to DOTcvpSB (dotcvp_sbml2dotcvpsb.m)
          9    Simulation
         10    GUI, a Graphical User Interface (dotcvp_gui.m)




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB   April 5, 2009   2 / 43
                                                                  Outline



   O UTLINE


          1    DOTcvpSB: Key Features
          2    Optimal Control Problems
          3    NLP and MINLP Solvers, Solution Scheme

          4    Single Optimization
          5    Multistart
          6    Sucessive Re-optimization
          7    Hybrid Strategy
          8    SBML to DOTcvpSB (dotcvp_sbml2dotcvpsb.m)
          9    Simulation
         10    GUI, a Graphical User Interface (dotcvp_gui.m)




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB   April 5, 2009   2 / 43
                                                           Brief information about the modules



   B RIEF INFORMATION ABOUT                                    THE MODULES

          1    Single Optimization – For most problems is this approach sufficient, but if the problem is multimodal, or
               high non-linear with many optimas, it is necessary to use an additional robust algorithm as, f.e. hybrid
               approach.
          2    Hybrid Strategy – Hybrid optimization is characterized by the combination of the deterministic local and
               the stochastic global method for the convergence speed up. Deterministic local methods stand out with
               their convergence speed into the nearest optimum, whereas the stochastic global methods are able to get
               to the vicinity of the global optimum without the extensive computational cost in comparison to deterministic
               global methods.
          3    Sucessive Re-optimization – Sucessive re-optimization is assigned to speed up the convergence for the
               problems with the high discretization level. At the beginning, the single optimization is run with relative low
               amount of the time intervals. It follows that the total number of optimized variables is lower. At the end of
               the optimization the amount of the time intervals is increased and the final values of the control variables
               are used for the next single optimization initialization. This iterative approach speeds up the final
               convergence of all optimizations.
          4    SBML to DOTcvpSB – This module is able to import SBML models defined as XML files.
          5    Simulation – Simulation module serves for problem simulation and state trajectories generation. It is useful
               to use this module in combination, f.e. with sbml2dotcvp module to input check.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                       Brief Tutorial on Using DOTcvpSB      April 5, 2009   3 / 43
                                                           Brief information about the modules



   B RIEF INFORMATION ABOUT                                    THE MODULES

          1    Single Optimization – For most problems is this approach sufficient, but if the problem is multimodal, or
               high non-linear with many optimas, it is necessary to use an additional robust algorithm as, f.e. hybrid
               approach.
          2    Hybrid Strategy – Hybrid optimization is characterized by the combination of the deterministic local and
               the stochastic global method for the convergence speed up. Deterministic local methods stand out with
               their convergence speed into the nearest optimum, whereas the stochastic global methods are able to get
               to the vicinity of the global optimum without the extensive computational cost in comparison to deterministic
               global methods.
          3    Sucessive Re-optimization – Sucessive re-optimization is assigned to speed up the convergence for the
               problems with the high discretization level. At the beginning, the single optimization is run with relative low
               amount of the time intervals. It follows that the total number of optimized variables is lower. At the end of
               the optimization the amount of the time intervals is increased and the final values of the control variables
               are used for the next single optimization initialization. This iterative approach speeds up the final
               convergence of all optimizations.
          4    SBML to DOTcvpSB – This module is able to import SBML models defined as XML files.
          5    Simulation – Simulation module serves for problem simulation and state trajectories generation. It is useful
               to use this module in combination, f.e. with sbml2dotcvp module to input check.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                       Brief Tutorial on Using DOTcvpSB      April 5, 2009   3 / 43
                                                           Brief information about the modules



   B RIEF INFORMATION ABOUT                                    THE MODULES

          1    Single Optimization – For most problems is this approach sufficient, but if the problem is multimodal, or
               high non-linear with many optimas, it is necessary to use an additional robust algorithm as, f.e. hybrid
               approach.
          2    Hybrid Strategy – Hybrid optimization is characterized by the combination of the deterministic local and
               the stochastic global method for the convergence speed up. Deterministic local methods stand out with
               their convergence speed into the nearest optimum, whereas the stochastic global methods are able to get
               to the vicinity of the global optimum without the extensive computational cost in comparison to deterministic
               global methods.
          3    Sucessive Re-optimization – Sucessive re-optimization is assigned to speed up the convergence for the
               problems with the high discretization level. At the beginning, the single optimization is run with relative low
               amount of the time intervals. It follows that the total number of optimized variables is lower. At the end of
               the optimization the amount of the time intervals is increased and the final values of the control variables
               are used for the next single optimization initialization. This iterative approach speeds up the final
               convergence of all optimizations.
          4    SBML to DOTcvpSB – This module is able to import SBML models defined as XML files.
          5    Simulation – Simulation module serves for problem simulation and state trajectories generation. It is useful
               to use this module in combination, f.e. with sbml2dotcvp module to input check.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                       Brief Tutorial on Using DOTcvpSB      April 5, 2009   3 / 43
                                                           Brief information about the modules



   B RIEF INFORMATION ABOUT                                    THE MODULES

          1    Single Optimization – For most problems is this approach sufficient, but if the problem is multimodal, or
               high non-linear with many optimas, it is necessary to use an additional robust algorithm as, f.e. hybrid
               approach.
          2    Hybrid Strategy – Hybrid optimization is characterized by the combination of the deterministic local and
               the stochastic global method for the convergence speed up. Deterministic local methods stand out with
               their convergence speed into the nearest optimum, whereas the stochastic global methods are able to get
               to the vicinity of the global optimum without the extensive computational cost in comparison to deterministic
               global methods.
          3    Sucessive Re-optimization – Sucessive re-optimization is assigned to speed up the convergence for the
               problems with the high discretization level. At the beginning, the single optimization is run with relative low
               amount of the time intervals. It follows that the total number of optimized variables is lower. At the end of
               the optimization the amount of the time intervals is increased and the final values of the control variables
               are used for the next single optimization initialization. This iterative approach speeds up the final
               convergence of all optimizations.
          4    SBML to DOTcvpSB – This module is able to import SBML models defined as XML files.
          5    Simulation – Simulation module serves for problem simulation and state trajectories generation. It is useful
               to use this module in combination, f.e. with sbml2dotcvp module to input check.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                       Brief Tutorial on Using DOTcvpSB      April 5, 2009   3 / 43
                                                           Brief information about the modules



   B RIEF INFORMATION ABOUT                                    THE MODULES

          1    Single Optimization – For most problems is this approach sufficient, but if the problem is multimodal, or
               high non-linear with many optimas, it is necessary to use an additional robust algorithm as, f.e. hybrid
               approach.
          2    Hybrid Strategy – Hybrid optimization is characterized by the combination of the deterministic local and
               the stochastic global method for the convergence speed up. Deterministic local methods stand out with
               their convergence speed into the nearest optimum, whereas the stochastic global methods are able to get
               to the vicinity of the global optimum without the extensive computational cost in comparison to deterministic
               global methods.
          3    Sucessive Re-optimization – Sucessive re-optimization is assigned to speed up the convergence for the
               problems with the high discretization level. At the beginning, the single optimization is run with relative low
               amount of the time intervals. It follows that the total number of optimized variables is lower. At the end of
               the optimization the amount of the time intervals is increased and the final values of the control variables
               are used for the next single optimization initialization. This iterative approach speeds up the final
               convergence of all optimizations.
          4    SBML to DOTcvpSB – This module is able to import SBML models defined as XML files.
          5    Simulation – Simulation module serves for problem simulation and state trajectories generation. It is useful
               to use this module in combination, f.e. with sbml2dotcvp module to input check.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                       Brief Tutorial on Using DOTcvpSB      April 5, 2009   3 / 43
                                                              DOTcvpSB       Toolbox Key Features



   T OOLBOX K EY F EATURES



               handling of a wide class of dynamic optimization problems, including constrained, unconstrained, fixed, and
               free terminal time problems described by ordinary differential equations (ODEs), as well as continuous and
               mixed integer decision variables;
               the inner initial value problem (IVP) is solved using the state-of-the-art methods available in SUNDIALS;
               the outer (MI)NLP problem can be solved using a number of advanced solvers, including local deterministic
               methods, stochastic global optimization methods, and hybrid metaheuristics;
               in addition to the traditional single optimization approach, the toolbox also offers more sophisticated
               strategies, like multistart, sucessive re-optimization, and hybrid strategies;
               a graphical user interface (GUI) which makes the definition and edition of a problem more easy and clear;
               possibility of importing SBML models;
               many output options for the results, including detailed figures.




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB                        April 5, 2009   4 / 43
                                                                 Optimal Control Problems    Problem Definition



   P ROBLEM D EFINITION
    Consider a dynamical system described by the vector of ordinary differential equations (ODEs)
                                                                              ˙
                                                                              x = f (t, x, u, p)
    with given initial and terminal conditions
                                                                    x(0) = x0 ,              x(tF ) = xF
    The aim is to find the optimal control policy ui , vector of the time-independent parameters p, and the final time tF ,
    when minimum or free time problem is considered that minimizes the cost function J0
                                                                   min J0 = G0 (tF , xF , u, p)
                                                                  tF ,ui ,p

    subject to the constraints, introduced also in general Mayer form as the cost function
                                                           Jl = Gl (tF , xF , u, p),              l = 1, me + mi
    where me and mi represent the equality and inequality constraints, respectively. The bounds on the optimized
    variables are given by
                                                              ∆ti ∈ [∆tilower bound , ∆tiupper bound ]
                                                                ui ∈ [uilower bound , uiupper bound ]
                                                                p ∈ [plower bound , pupper bound ]
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                   Brief Tutorial on Using DOTcvpSB        April 5, 2009   5 / 43
                                                                 NLP and MINLP solvers     MI/NLP Formulation



   S OLUTION OF THE R ESULTING NLP OR MINLP P ROBLEMS

    The dynamic optimization (DO) problem can be translated into a set of optimized variables

                                                                                       T
                                                           y T = [∆t1 , . . . , ∆tN , u1 , . . . , uN , p]
                                                                                                    T



    that can be handled by a suitable NLP or Mixed-Integer (MI) NLP solver.

    The toolbox is able to interface following MI/NLP solvers
         deterministic
       IPOPT: Interior Point OPTimizer
    FMINCON: Find MINimum of CONstrained nonlinear multivariable function
      MISQP: Mixed-Integer Sequential Quadratic Programming
               stochastic
              DE: Differential Evolution
            SRES: Stochastic Ranking Evolution Strategy
               and hybrid
          ACO MI : Ant Colony Optimization for Mixed Integer nonlinear programming problems
           MITS: Mixed-Integer Tabu Search algorithm



Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()                 Brief Tutorial on Using DOTcvpSB       April 5, 2009   6 / 43
                                                                   DOTcvpSB       Scheme of the Solution of the Resulting NLP Problem



   A LGORITHM S CHEME




                F IGURE : DOTcvpSB solution scheme of the resulting dynamic optimization problem with sensitivity equations approach.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()        Brief Tutorial on Using DOTcvpSB                                        April 5, 2009   7 / 43
                                                                         DOTcvpSB       Organization of the Toolbox Code



   C ODE O RGANIZATION




                                                           F IGURE : Organization of the toolbox code.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()              Brief Tutorial on Using DOTcvpSB                     April 5, 2009   8 / 43
                                                           Single Optimization   Single Optimization



   S INGLE O PTIMIZATION

      Single Optimization



    T HIS PROBLEM CAN BE RUN FROM THE TOOLBOX WITH THE FOLLOWING NAME
    cdop_DrugDisplacementProblemA.m                     < input file with all options
    midop_PhaseResettingOfCalciumOscillationsA.m
    cdop_DrugDisplacementProblemA_simple.m              < input file with minimum
    midop_PhaseResettingOfCalciumOscillationsA_simple.m   options, undefined options
                                                          are taken from the default
                                                          file
    gui_DrugDisplacementProblemA.dotcvp                 < input file defined by the
                                                          DOTcvp graphical user
                                                          interface (GUI)

    T HE FILE WITH THE DEFAULT SETTINGS FOR THE PRESENTED MODULE IS THE FOLLOWING
    dotcvp_single_optimization_default.m


Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()       Brief Tutorial on Using DOTcvpSB        April 5, 2009   9 / 43
                                                           Single Optimization   Single Optimization



   S INGLE O PTIMIZATION

      Single Optimization



    T HIS PROBLEM CAN BE RUN FROM THE TOOLBOX WITH THE FOLLOWING NAME
    cdop_DrugDisplacementProblemA.m                     < input file with all options
    midop_PhaseResettingOfCalciumOscillationsA.m
    cdop_DrugDisplacementProblemA_simple.m              < input file with minimum
    midop_PhaseResettingOfCalciumOscillationsA_simple.m   options, undefined options
                                                          are taken from the default
                                                          file
    gui_DrugDisplacementProblemA.dotcvp                 < input file defined by the
                                                          DOTcvp graphical user
                                                          interface (GUI)

    T HE FILE WITH THE DEFAULT SETTINGS FOR THE PRESENTED MODULE IS THE FOLLOWING
    dotcvp_single_optimization_default.m


Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()       Brief Tutorial on Using DOTcvpSB        April 5, 2009   9 / 43
                                                                   Single Optimization   Drug Displacement Problem without the path constraint



   E XAMPLE 1: D RUG D ISPLACEMENT P ROBLEM WITHOUT THE PATH CONSTRAINT




                                              Example 1: Drug Displacement Problem without the path constraint




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga        ()          Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   10 / 43
                                                                        Single Optimization    Drug Displacement Problem without the path constraint



   P ROBLEM D EFINITION
    The problem consists of the right rate projection of phenylbutazone infusion to minimize the time needed to reach
    in a patient’s bloodstream a desired level of two drugs
                                                                                  min J0 {tF }
                                                                                   ui ,ti

    subject to
                                                            ˙
                                                            x1 = g4 (g3 (0.02 − x1 ) + 46.4x1 (u − 2x2 ))
                                                            ˙ 2 = g4 (g2 (u − 2x2 )) + 46.4(0.02 − x1 )
                                                            x

    with gi , i = 1, 4 defined as follows
                                                                                                                                                                   2
                                                            2                                       2                                                             g1
          g1 = 1 + 0.2(x1 + x2 ),                     g2 = g1 + 232 + 46.4x2 ,                g3 = g1 + 232 + 46.4x1 ,                            g4 =
                                                                                                                                                         g2 g3 − 2152.96x1 x2
    where the state variables represent the concentration of warfarin and phenylbutazone drugs. The initial values
    were set for the process and decisions variables at the value of x0 = [0.02; 2.00] and u0 = [4], respectively. The
    boundaries of the decision variables are as follows: u ∈ [0; 8].
    The equality point constraints on the final amount of displacement drugs are given by
                                                                   x1 (tF ) = 0.02,             x2 (tF ) = 2.00

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga           ()            Brief Tutorial on Using DOTcvpSB                                                          April 5, 2009   11 / 43
                                                                        Single Optimization   Drug Displacement Problem without the path constraint



   U SER I NPUT: P ROBLEM AND IVP I NITIALIZATION

         1   clear mex; clear all; close all; global data;
         2   % --------------------------------------------------- %
         3   % Initialization:
         4   % --------------------------------------------------- %
         5   data.name                 = ’DrugDisplacementProblemA_simple’;
         6
         7   % --------------------------------------------------- %
         8   % Settings for IVP (ODEs, sensitivities)
         9   % --------------------------------------------------- %
        10   data.odes.res(1)          = {’((1+0.2*(y(1)+y(2)))^2/(((1+0.2*(y(1)+y(2)))^2+232+46.4*y(2)) *((1+0.2*(y(1)+y(2)))^2+232+46.4*y(1))
        11    -2152.96*y(1)*y(2)))*(((1+0.2*(y(1)+y(2)))^2+232+46.4*y(1))*(0.02-y(1))+46.4*y(1)*(u(1)-2*y(2)))’};
        12   data.odes.res(2)          = {’((1+0.2*(y(1)+y(2)))^2/(((1+0.2*(y(1)+y(2)))^2+232+46.4*y(2)) *((1+0.2*(y(1)+y(2)))^2+232+46.4*y(1))
        13    -2152.96*y(1)*y(2)))*(((1+0.2*(y(1)+y(2)))^2+232+46.4*y(2))*(u(1)-2*y(2))+46.4*(0.02-y(1)))’};
        14   data.odes.res(3)          = {’1’};
        15   data.odes.ic              = [0.02 0.0 0.0];
        16   data.odes.tf              = 500.0; %final time
        17   data.odes.RelTol          = 1*10^(-8); %IVP relative tolerance level
        18   data.odes.AbsTol          = 1*10^(-8); %IVP absolute tolerance level
        19   data.sens.SensAbsTol      = 1*10^(-8); %absolute tolerance for sensitivity variables



               DATA . NAME         The variable which contains the name of the problem.
 DATA . ODES . RES (1–3)           Here is defined the mathematical description of the system in the standard ODEs form. The notation y(number),
                                   u(number), p(number) has to be used for the state, decision variables, and time-independent parameters,
                                   respectively.
             DATA . ODES . IC      The vector of the initial state variables defined for the system.
             DATA . ODES . TF      The final time of the optimization, initial time is set as default at the value of 0.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga         ()              Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   12 / 43
                                                                       Single Optimization   Drug Displacement Problem without the path constraint



   U SER I NPUT: NLP D EFINITION


         1    % --------------------------------------------------- %
         2    % NLP definition:
         3    % --------------------------------------------------- %
         4    data.nlp.RHO              = 5; %number of time intervals
         5    data.nlp.J0               = ’y(3)’; %cost function: min-max(cost function)
         6    data.nlp.u0               = [4.0]; %initial value for control values
         7    data.nlp.lb               = [0.0]; %lower bounds for control values
         8    data.nlp.ub               = [8.0]; %upper bounds for control values
         9    data.nlp.solver           = ’IPOPT’; %[’FMINCON’|’IPOPT’|’SRES’|’DE’|’ACOMI’|’MISQP’|’MITS’]
        10    data.nlp.FreeTime         = ’on’; %[’on’|’off’] set ’on’ if free time is considered



             DATA . NLP.RHO        The number of time intervals.
               DATA . NLP.J0       The cost function is defined in the standard Mayer form. The problem is set as default to find a minimum of the
                                   presented cost function.
               DATA . NLP. U 0     The initial values of the decision variables, if here is defined only one number for each control trajectory, then the
                                   initial control trajectory is constant, otherwise is needed to define, e.g. [u01 ; u02 ; ..] whether only one initial
                                   non-constant control trajectory is considered.
               DATA . NLP. LB      The vector of lower bounds of the decision variables.
               DATA . NLP. UB      The vector of upper bounds of the decision variables.



Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga        ()              Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   13 / 43
                                                                         Single Optimization   Drug Displacement Problem without the path constraint



   U SER I NPUT: E QUALITY C ONSTRAINTS

         1   % --------------------------------------------------- %
         2   % Equality constraints (ECs):
         3   % --------------------------------------------------- %
         4   data.nlp.eq.status        = ’on’; %[’on’|’off’] ECs
         5   data.nlp.eq.NEC           = 2; %number of active ECs
         6   data.nlp.eq.eq(1)         = {’y(1)-0.02’};
         7   data.nlp.eq.eq(2)         = {’y(2)-2.0’};
         8   data.nlp.eq.time(1)       = data.nlp.RHO;
         9   data.nlp.eq.time(2)       = data.nlp.RHO;



      DATA . NLP. EQ . STATUS          The information about the equality constrains, if those are active, set this option at the value of ’on’ otherwise ’off’.
         DATA . NLP. EQ .NEC           The number of active equality constraints.
     DATA . NLP. EQ . EQ (1-2)         The equality constraints defined by the user.
 DATA . NLP. EQ . TIME (1-2)           The segments after which are before defined equality constrains active.
                            N OTE The equality constraint of the type ’data.nl.eq.eq(1) = {’y(2)+0.1’}’ is active, the status is set at the value of ’on’ what
                                  means that the equation

                                                                                                     x2 (tF ) = −0.1

                                       is active in the final segment ’data.nlp.eq.time(1) = data.nlp.RHO’. This direction – timing is equal to the final time.
                                       If we would like to set the active equality constraint on the half time of the optimization, we have to define the before
                                       mentioned option on the value of 15 or ’data.nlp.RHO/2’.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga          ()              Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   14 / 43
                                                                           Single Optimization   Drug Displacement Problem without the path constraint



   U SER I NPUT: F INAL O PTIONS


         1   % --------------------------------------------------- %
         2   % Options for setting of the final output:
         3   % --------------------------------------------------- %
         4   data.options.intermediate = ’off’; %[’on’|’off’|’silent’] display of the intermediate results
         5   data.options.trajectories = 2; %how many state trajectories will be displayed
         6
         7   % --------------------------------------------------- %
         8   % Call of the main function (you do not change this!):
         9   % --------------------------------------------------- %
        10   dotcvp_main(data)



 DATA . OPTIONS . INTERMEDIATE                   The option related to the display of the intermediate results.
  DATA . OPTIONS . TRAJECTORIES                  How many state trajectories will be displayed in the final figure.
             DOT CVP _M AIN ( DATA ) Here the main function of the toolbox is called together with the data input structure.
                                      N OTE Only this file is needed for introducing a new problem into toolbox, everything else is generated
                                            automatically. Remember that this is only a simple input file. The undefined options are taken from the
                                            default settings [dotcvp_single_optimization_default.m], of course it is possible to change them. After
                                            MATLAB is started, it is necessary to run the [dotcvp_install.m] file to save the paths into MATLAB
                                            environment. Note, if you quit MATLAB environment, all paths will be deleted, so you need to run the
                                            [dotcvp_install.m] file with the new start of MATLAB again.



Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga             ()             Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   15 / 43
                                                                 Single Optimization     Drug Displacement Problem without the path constraint



   S CREEN O UTPUT: I NITIALIZATION

         1   ________________________________________________________
         2
         3    DOTcvp: Dynamic Optimization Toolbox with CVP approach
         4     for handling continuous and mixed-integer DO problems
         5
         6     Main author: Tomas Hirmajer, thirmajer@iim.csic.es
         7       Coauthors: Eva Balsa-Canto and Julio R. Banga
         8       Web pages: http://www.iim.csic.es/~dotcvp/
         9                  http://www.iim.csic.es/~dotcvpsb/
        10    Core version: DOTcvp_R2008_E2
        11   ________________________________________________________
        12
        13   ________________________________________________________
        14
        15    DOTcvp - a Module for Single Optimization
        16   ________________________________________________________
        17
        18   Saving of the ODE ................................................. done!
        19   Generation of the file: cvm_rhs.m (ODE - MATLAB) .................. done!
        20   Saving of the parameters .......................................... done!
        21   Generation of the file: cvm_d/bjac.m (Jacobian - MATLAB) .......... done!
        22   Saving of the cost function (J0) .................................. done!
        23   Generation of the file: temp_TimeDiscontinuity .................... done!
        24   Generation of the gradients (J0) .................................. done!
        25   Generation of the file: cvm_rhsS.m (sensitivities - MATLAB) ....... done!
        26   Saving of the in/equality constraints (Ji) ........................ done!
        27   Generation of the gradients (Ji) .................................. done!
        28   Generation of the file: temp_cvfdx.m (main IVP file) .............. done!
        29   Optimizing of the process (N=5; min(J0); IPOPT; DrugDisplacementProblemA_simple) ... in progress
        30   Default settings are loading ...................................... done!


Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()               Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   16 / 43
                                                                     Single Optimization   Drug Displacement Problem without the path constraint



   S CREEN O UTPUT: R ESULTS

         1   ____________________________
         2   Final results [single-optimization]:
         3    ............. Problem name: DrugDisplacementProblemA_simple
         4    ...... NLP or MINLP solver: IPOPT
         5    . Number of time intervals: 5
         6    ... IVP relative tolerance: 1.000000e-008
         7    ... IVP absolute tolerance: 1.000000e-008
         8    . Sens. absolute tolerance: 1.000000e-008
         9    ............ NLP tolerance: 1.000000e-005
        10    ....... Final state values: 2.000004e-002 2.000014e+000 2.212422e+002
        11    ...... 1th optimal control: 7.999963e+000 7.999963e+000 7.999963e+000 8.000000e+000 5.454379e-007
        12    ..... Final size of the dt: 1.000009e-001 1.000009e-001 1.000009e-001 1.886241e+002 3.231812e+001
        13    ..... Final time [sum(dt)]: 2.212422e+002
        14   ____________________________
        15    ............ Final CPUtime: 18.89062500 seconds
        16    . Cost function [min(J_0)]: 221.24223300
        17
        18    The detailed information is saved to the workspace structure with the name ’data’.
        19
        20   data =
        21            name:      ’DrugDisplacementProblemA_simple’
        22            odes:      [1x1 struct]
        23            sens:      [1x1 struct]
        24             nlp:      [1x1 struct]
        25         options:      [1x1 struct]
        26         version:      ’DOTcvp_R2008_E2’
        27        compiler:      ’None’
        28          output:      [1x1 struct]
        29               p:      [5.4544e-007 5]
        30          gradJ0:      [0 0 0 0 0 1 1 1 1 1]
        31          gradJi:      [2x10 double]
        32              J0:      221.2422
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga         ()           Brief Tutorial on Using DOTcvpSB                                          April 5, 2009   17 / 43
                                                                                Single Optimization                    Drug Displacement Problem without the path constraint



   G RAPHICAL O UTPUT: O PTIMAL T RAJECTORIES

                        0.03                                                                          4
                                   x1
                                                                                                                                              8
                                   x2


                                                                                                                                              7


                                                                                                                                              6


                                                                                                                                              5




                                                                                                                           Control variable
      State Variable




                                                                                                          State Variable
                       0.025
                                                                                                      2                                       4


                                                                                                                                              3


                                                                                                                                              2


                                                                                                                                              1

                        0.02
                                                                                                                                              0
                                                                                                      0
                               0          50           100           150           200                                                            0   50             100          150   200
                                                             Time                                                                                                          Time



                                   F IGURE : Optimal state trajectories (left) and the optimal control profile (right) for the drug displacement problem.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                 ()              Brief Tutorial on Using DOTcvpSB                                                                    April 5, 2009   18 / 43
                                                                Single Optimization   Resetting of Calcium Oscillator



   E XAMPLE 2: R ESETTING OF C ALCIUM O SCILLATOR




                                                      Example 2: Resetting of Calcium Oscillator




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga     ()          Brief Tutorial on Using DOTcvpSB                    April 5, 2009   19 / 43
                                                                 Single Optimization   Resetting of Calcium Oscillator



   P ROBLEM D EFINITION
    The model describes intracellular calcium spiking in hepatocytes induced by an extracellular increase in
    adenosine triphosphate concentration. The problem is formulated as minimization of the state variables
    deviations (oscillatory behaviour) from the desired values
                                                                                           
                                          22    4                                               
                                                                 s 2
                                min J0            wj xj (t) − xj    + w5 u1 (t) + w6 u2 (t) dt
                             x(t),u(t),p  0                                                     
                                                           j=1

    subject to
                                                     k3 x1 x2     k5 x1 x3
                                  ˙
            activated G-protein → x1 = k1 + k2 x1 −            −
                                                    x1 + K4      x1 + K6
                                                           k8 x2
                                ˙
    active phospholipase C → x2 = (1 − u2 )k7 x1 −
                                                        x2 + K9
                                     k10 x2 x3 x4                          k16 x3   x4         k14 x3                 k14 x3
                                ˙
        intracellular calcium → x3 =              + k12 x2 + k13 x1 −             +    − u1             − (1 − u1 )
                                     x4 + K11                            x3 + K17   10      p1 x3 + K15             x3 + K15
                                       k10 x2 x3 x4      k16 x3       x4
                                ˙
            intra-ER calcium → x4 = −               +             −
                                       x4 + K11        x3 + K17       10
    and the time-independent parameter: 1 ≤ p1 ≤ 1.3
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()             Brief Tutorial on Using DOTcvpSB                    April 5, 2009   20 / 43
                                                                                      Single Optimization   Resetting of Calcium Oscillator



   S IMULATION WITH NO AND F ULL I NHIBITION

                          30                                                                                                   20
                                x                                                                                                    x
                                    1                                                                                                    1
                                x                                                                                              18    x
                                    2                                                                                                    2

                          25    x                                                                                                    x
                                    3                                                                                                    3
                                x                                                                                              16    x
                                    4                                                                                                    4

                                                                                                                               14
                          20
                                                                                                                               12
        State variables




                                                                                                             State variables
                          15                                                                                                   10

                                                                                                                                8
                          10
                                                                                                                                6

                                                                                                                                4
                           5
                                                                                                                                2

                           0                                                                                                    0
                            0   2       4   6   8     10     12   14        16   18       20       22                            0   2       4   6   8   10     12   14   16   18      20       22
                                                        Time                                                                                               Time


    F IGURE : Simulation of the system with no inhibition (left: u1 = 0, u2 = 0, p1 = 1) and with the constant maximum inhibition of the
    PMCA (plasma membrane Ca2+ ) ion pump (right: u1 = 1, u2 = 0, p1 = 1.3) for the calcium oscillator problem.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                    ()                 Brief Tutorial on Using DOTcvpSB                                                          April 5, 2009    21 / 43
                                                                                      Single Optimization   Resetting of Calcium Oscillator



   T HE S CENARIO WITH F REE T RANSITION T IMES AND ONE C ONTROL VARIABLE

                          25
                                 x                                                                                              1
                                     1
                                 x
                                     2                                                                                         0.9
                                 x
                                     3
                          20     x                                                                                             0.8
                                     4


                                                                                                                               0.7




                                                                                                            Control variable
                          15                                                                                                   0.6
        State variables




                                                                                                                               0.5


                          10                                                                                                   0.4

                                                                                                                               0.3

                                                                                                                               0.2
                           5
                                                                                                                               0.1

                                                                                                                                0
                           0
                            0    2       4   6   8    10     12   14        16   18       20       22                            0   2   4    6   8   10     12   14   16   18       20          22
                                                        Time                                                                                            Time


                                F IGURE : Optimal state trajectories (left) with corresponding control profile (right) for the calcium oscillator problem.


Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                    ()                 Brief Tutorial on Using DOTcvpSB                                                       April 5, 2009        22 / 43
                                                                                                               Single Optimization                            Resetting of Calcium Oscillator



   T HE S CENARIO WITH F REE T RANSITION T IMES AND T WO C ONTROL VARIABLES

                                                                                       Activated G−proteins                                         Active Phospholipase C                                  Pump Inhibitor
                                                                          12                                                           30
                                                                                                                                                                                                   1
                                                                          10                                                           25
                                                                                                                                                                                                  0.8




                                                                                                                                                                               Control variable
                                                      State variable




                                                                                                                   State variable
                                                                               8                                                       20
                                                                                                                                                                                                  0.6
                                                                               6                                                       15
                                                                                                                                                                                                  0.4
                                                                               4                                                       10
                                                                                                                                                                                                  0.2
                                                                               2                                                            5
                                                                                                                                                                                                   0
                                                                               0                                                            0
                                                                                   0            10         20                                   0           10           20                             0        10           20
                                                                                                Time                                                        Time                                                 Time

                                                                                       Intracellular Calcium                                          Intra−ER Calcium                                      Channel Blocker
                                                                               5                                                            4
                                                                                                                                                                                                   1

                                                                               4                                                                                                                  0.8
                                                                                                                                            3




                                                                                                                                                                               Control variable
                                                              State variable




                                                                                                                           State variable
                                                                               3                                                                                                                  0.6
                                                                                                                                            2
                                                                               2                                                                                                                  0.4

                                                                                                                                            1                                                     0.2
                                                                               1

                                                                                                                                                                                                   0
                                                                               0                                                            0
                                                                                   0            10         20                                   0           10           20                             0        10           20
                                                                                                Time                                                        Time                                                 Time



                                              F IGURE : Optimal state and control trajectories for the calcium oscillator problem.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                                        ()                      Brief Tutorial on Using DOTcvpSB                                                                                April 5, 2009   23 / 43
                                                                                           Single Optimization     Resetting of Calcium Oscillator



   T HE P ROBLEM S IMULATION TO L ONGER T IMES
                                                                        30
                                                                                                                                                          x1
                                                                                                                                                          x2
                                                                                                                                                          x3
                                                                                                                                                          x4
                                                                        25




                                                                        20




                                                      State variables   15




                                                                        10




                                                                         5




                                                                         0
                                                                             0        10    20              30             40             50         60
                                                                                                          Time


                     F IGURE : The problem simulation to longer times (tF = 60) without more stimuli for the calcium oscillator problem.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                              ()              Brief Tutorial on Using DOTcvpSB                              April 5, 2009   24 / 43
                                                           Hybrid Strategy   Hybrid Strategy



   H YBRID S TRATEGY


      Hybrid Strategy


    W HAT IS NEEDED TO BE DONE IS JUST TO ADD OR CHANGE THE OPTION IN THE USER INPUT FILE AS FOLLOWS
    data.options.action=’hybrid-strategy’

    T HIS PROBLEM CAN BE RUN FROM THE TOOLBOX WITH THE FOLLOWING NAME
    cdop_DrugDisplacementProblemB.m        < input file with all options
    cdop_DrugDisplacementProblemB_simple.m < input file with minimum options, undefined
                                            options are taken from the default file

    T HE FILE WITH THE DEFAULT SETTINGS FOR THE PRESENTED MODULE IS THE FOLLOWING
    dotcvp_hybrid_strategy_default.m




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB     April 5, 2009   25 / 43
                                                           Hybrid Strategy   Hybrid Strategy



   H YBRID S TRATEGY


      Hybrid Strategy


    W HAT IS NEEDED TO BE DONE IS JUST TO ADD OR CHANGE THE OPTION IN THE USER INPUT FILE AS FOLLOWS
    data.options.action=’hybrid-strategy’

    T HIS PROBLEM CAN BE RUN FROM THE TOOLBOX WITH THE FOLLOWING NAME
    cdop_DrugDisplacementProblemB.m        < input file with all options
    cdop_DrugDisplacementProblemB_simple.m < input file with minimum options, undefined
                                            options are taken from the default file

    T HE FILE WITH THE DEFAULT SETTINGS FOR THE PRESENTED MODULE IS THE FOLLOWING
    dotcvp_hybrid_strategy_default.m




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB     April 5, 2009   25 / 43
                                                           Hybrid Strategy   Hybrid Strategy



   H YBRID S TRATEGY


      Hybrid Strategy


    W HAT IS NEEDED TO BE DONE IS JUST TO ADD OR CHANGE THE OPTION IN THE USER INPUT FILE AS FOLLOWS
    data.options.action=’hybrid-strategy’

    T HIS PROBLEM CAN BE RUN FROM THE TOOLBOX WITH THE FOLLOWING NAME
    cdop_DrugDisplacementProblemB.m        < input file with all options
    cdop_DrugDisplacementProblemB_simple.m < input file with minimum options, undefined
                                            options are taken from the default file

    T HE FILE WITH THE DEFAULT SETTINGS FOR THE PRESENTED MODULE IS THE FOLLOWING
    dotcvp_hybrid_strategy_default.m




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()   Brief Tutorial on Using DOTcvpSB     April 5, 2009   25 / 43
                                                                          Hybrid Strategy    Hybrid Strategy



   S ETTINGS F OR T HE H YBRID S TRATEGY

         1   function [data] = dotcvp_hybrid_strategy_default(data)
         2
         3   switch data.option.Hybrid_method
         4
         5         case{’stochastic’}
         6             % first step [a stochastic solver]: please fill basic settings for
         7             % a stochastic solver, otherwise they will be taken from the user
         8             % input file
         9
        10              data.option.HybridStrategy.IVPRelTol      =   1e-005;   %   IVP relative tolerance level
        11              data.option.HybridStrategy.IVPAbsTol      =   1e-005;   %   IVP absolute tolerance level
        12              data.option.HybridStrategy.NLPTol         =   1e-003;   %   NLP tolerance level
        13              data.option.HybridStrategy.NLPSolver      =   ’DE’;     %   chose a solver [’FMINCON’|’IPOPT’|’SRES’|’DE’|’ACOMI’|’MISQP’|’MITS’]
        14              data.option.HybridStrategy.NLPsettings    =   ’None’;   %   insert the file name that contains settings for NLP solver, otherwise ’None’
        15              data.option.HybridStrategy.MaxIter        =   50;       %   maximum number of iterations
        16              data.option.HybridStrategy.MaxCPUTime     =   inf;      %   maximum CPU time of the optimization (60*60*0.25) = 15 minutes
        17              data.option.HybridStrategy.intermediate   =   ’off’;    %   [’on’|’off’|’silent’] display of the intermediate results
        18
        19         case{’deterministic’}
        20             % second step [a deterministic solver]: all settings for a
        21             % deterministic solver are taken from the user input file, but the
        22             % user can introduce some new or special settings, which rewrites
        23             % those from the input file. The name of the structure is the same
        24             % as the name in the user input file, e.g. data.nlp.RHO, etc.
        25
        26   end
        27
        28   end




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga      ()                Brief Tutorial on Using DOTcvpSB                                         April 5, 2009   26 / 43
                                                                       Hybrid Strategy   Drug Displacement Problem with the path constraint



   E XAMPLE 3: D RUG D ISPLACEMENT P ROBLEM WITH THE PATH CONSTRAINT




                                                 Example 3: Drug Displacement Problem with the path constraint




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga          ()        Brief Tutorial on Using DOTcvpSB                                       April 5, 2009   27 / 43
                                                                            Hybrid Strategy    Drug Displacement Problem with the path constraint



   P ROBLEM D EFINITION
    The problem consists of the right rate projection of phenylbutazone infusion to minimize the time needed to reach
    in a patient’s bloodstream a desired level of two drugs
                                                                                  min J0 {tF }
                                                                                   ui ,ti

    subject to
                                                            ˙
                                                            x1 = g1 (g4 (0.02 − x1 ) + 46.4x1 (u − 2x2 ))
                                                            ˙ 2 = g1 (g3 (u − 2x2 )) + 46.4(0.02 − x1 )
                                                            x

    with gi , i = 1, 4 defined as follows
                                                                                                                                                                     2
                                                            2                                       2                                                               g2
          g2 = 1 + 0.2(x1 + x2 ),                     g3 = g2 + 232 + 46.4x2 ,                g4 = g2 + 232 + 46.4x1 ,                              g1 =
                                                                                                                                                           g3 g4 − 2152.96x1 x2
    where the state variables represent the concentration of warfarin and phenylbutazone drugs. The initial values
    were set for the process and decisions variables at the value of x0 = [0.02; 2.00] and u0 = [4], respectively. The
    boundaries of the decision variables are as follows: u ∈ [0; 8].
    The equality point constraints on the final amount of displacement drugs are given by
                                                         x1 (tF ) = 0.02,      x2 (tF ) = 2.00,                 x1 (t) ≤ 0.026

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga           ()            Brief Tutorial on Using DOTcvpSB                                                            April 5, 2009   28 / 43
                                                                              Hybrid Strategy    Drug Displacement Problem with the path constraint



   U SER I NPUT: I NEQUALITY C ONSTRAINTS AND SOME F INAL O; PTIONS

         1   % --------------------------------------------------- %
         2   % Inequality /path/ constraints (INECs):
         3   % --------------------------------------------------- %
         4   data.nlp.ineq.status      = ’on’; %[’on’|’off’] INECs
         5   data.nlp.ineq.NEC         = 1; %number of active INECs
         6   data.nlp.ineq.InNUM       = 0; %how many inequality constraints are ’>’ else ’<’
         7   data.nlp.ineq.eq(1)       = {’y(1)-0.026’};
         8   data.nlp.ineq.Tol         = 10^(-6); %tolerance level of violation of INECs
         9   data.nlp.ineq.PenaltyFun = ’on’; %[’on’|’off’] INECs penalty function
        10   data.nlp.ineq.PenaltyCoe = 10^(5); %J0=J0+data.nlp.ineq.PenaltyCoe*ViolationOfInequalityConstraint
        11                                        /* for every inequality constraint one parameter */
        12
        13   % --------------------------------------------------- %
        14   % Options for setting of the final output:
        15   % --------------------------------------------------- %
        16   data.options.trajectories = 2; %how many state trajectories will be displayed
        17
        18   data.options.action                  = ’hybrid-strategy’; %[’single-optimization’|’re-optimization’|’hybrid-strategy’|’simulation’]



  DATA . NLP. INEQ . STATUS            The information about the inequality constrains, if those are active, set this option at the value of ’on’ otherwise ’off’.
      DATA . NLP. INEQ .NEC            The number of active inequality constraints.
 DATA . NLP. INEQ .I N NUM             In this option the user sets how many inequality constraints are ’>’ else ’<’.
     DATA . NLP. INEQ . EQ (1)         The inequality constraints defined by the user.
       DATA . NLP. INEQ .T OL          The tolerance level of violation of the inequality constraints.
   DATA . OPTIONS . ACTION             The hybrid strategy was chosen instead of a single optimization module.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga              ()            Brief Tutorial on Using DOTcvpSB                                       April 5, 2009   29 / 43
                                                                            Hybrid Strategy             Drug Displacement Problem with the path constraint



   G RAPHICAL O UTPUT: O PTIMAL T RAJECTORIES

                                                                                 4
                                       x1                                                                                        8
                                       x2
                                                                                 3.5
                                                                                                                                 7


                                                                                 3                                               6
                           0.027


                           0.026                                                 2.5                                             5




                                                                                                              Control variable
          State Variable




                                                                                       State Variable
                           0.025
                                                                                 2                                               4
                           0.024
                                                                                                                                 3
                           0.023                                                 1.5


                           0.022                                                                                                 2
                                                                                 1

                           0.021                                                                                                 1
                                                                                 0.5
                            0.02
                                                                                                                                 0

                           0.019                                                 0
                                   0        50   100      150   200   250                                                            0   50     100             150   200   250
                                                       Time                                                                                                  Time



    F IGURE : Optimal state trajectories (left) and the optimal control profile (right) for the drug displacement problem with the path constraint.

    Final value of the cost function is 266.0897 (tF ) and violation of the constraints is less than 10−8 .
    The Stochastic part is secured by DE and the deterministic part by MISQP.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga             ()          Brief Tutorial on Using DOTcvpSB                                                                 April 5, 2009   30 / 43
                                                           SBML to DOTcvpSB      SBML to DOTcvpSB



   I MPORTING SBML M ODELS




                                                           Importing SBML Models




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()       Brief Tutorial on Using DOTcvpSB     April 5, 2009   31 / 43
                                                              SBML to DOTcvpSB      SBML to DOTcvpSB



   SBML TO DOT CVP SB


      SBML to DOTcvpSB



    T HE TRANSLATOR FROM THE SBML,                    TYPICALLY THE ’ XML’ FILES TO THE                DOT CVP SB ’. DOTCVP ’ CAN BE DONE BY
    TYPING THE FOLLOWING NAME
    dotcvp_sbml2dotcvpsb.m



    The problem can be downloaded as ’.xml’ file from the following web page:

    http://www.ebi.ac.uk/biomodels/

      BIOMD0000000005 Tyson1991_CellCycle_6var 1831270



Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()          Brief Tutorial on Using DOTcvpSB                                 April 5, 2009   32 / 43
                                                              SBML to DOTcvpSB      SBML to DOTcvpSB



   SBML TO DOT CVP SB


      SBML to DOTcvpSB



    T HE TRANSLATOR FROM THE SBML,                    TYPICALLY THE ’ XML’ FILES TO THE                DOT CVP SB ’. DOTCVP ’ CAN BE DONE BY
    TYPING THE FOLLOWING NAME
    dotcvp_sbml2dotcvpsb.m



    The problem can be downloaded as ’.xml’ file from the following web page:

    http://www.ebi.ac.uk/biomodels/

      BIOMD0000000005 Tyson1991_CellCycle_6var 1831270



Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()          Brief Tutorial on Using DOTcvpSB                                 April 5, 2009   32 / 43
                                                               SBML to DOTcvpSB      cdc2-cyclin Model



   E XAMPLE 4:                       CDC 2- CYCLIN         M ODEL




                                                            Example 4: cdc2-cyclin Model




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()           Brief Tutorial on Using DOTcvpSB      April 5, 2009   33 / 43
                                                                      SBML to DOTcvpSB      SBML to DOTcvpSB: Part 1



   SBML TO DOT CVP SB: PART 1

         1   ______________________________________________________________
         2
         3    DOTcvpSB: Dynamic Optimization Toolbox with CVP approach for
         4        handling continuous and MIDO problems in Systems Biology
         5
         6     Main author: Tomas Hirmajer, thirmajer@iim.csic.es
         7       Coauthors: Eva Balsa-Canto and Julio R. Banga
         8        Web page: http://www.iim.csic.es/~dotcvpsb/
         9    Core version: DOTcvp_R2008_E2
        10   _______________________________________________________________
        11
        12   _______________________________________________________________
        13
        14    DOTcvpSB - a Module for Import of SBML Models (.xml)
        15    Note: this module requires the installation of SBML toolbox.
        16   _______________________________________________________________
        17
        18   Do you want to print all the import information on the screen? [y|n]: y
        19
        20           Name:       Tyson1991_CellCycle_6var
        21             ID:       Tyson1991CellModel_6
        22     SBML_level:       2
        23   SBML_version:       1
        24    Time_symbol:
        25
        26   --- Getting information from model       functions ------------------
        27       1. parameter from the model          : k6 (k6_Reaction1) = 1.000000 (1.000000)
        28       2. parameter from the model          : k8notP (k8notP_Reaction2) = 1000000.000000 (1000000.000000)
        29       3. parameter from the model          : k9 (k9_Reaction3) = 1000.000000 (1000.000000)
        30       4. parameter from the model          : k3 (k3_Reaction4) = 200.000000 (200.000000)
        31       5. parameter from the model          : k5notP (k5notP_Reaction5) = 0.000000 (0.000000)
        32       6. parameter from the model          : k1aa (k1aa_Reaction6) = 0.015000 (0.015000)
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga         ()            Brief Tutorial on Using DOTcvpSB             April 5, 2009   34 / 43
                                                                            SBML to DOTcvpSB      SBML to DOTcvpSB: Part 2



   SBML TO DOT CVP SB: PART 2

         1        7.   parameter     from   the   model   :   k2 (k2_Reaction7) = 0.000000 (0.000000)
         2        8.   parameter     from   the   model   :   k7 (k7_Reaction8) = 0.600000 (0.600000)
         3        9.   parameter     from   the   model   :   k4 (k4_Reaction9) = 180.000000 (180.000000)
         4       10.   parameter     from   the   model   :   k4prime (k4prime_Reaction9) = 0.018000 (0.018000)
         5
         6   --- Getting information from reaction functions ---------------
         7       1. parameter for reaction( 1) : k6 (k6_Reaction1) = 1.000000 (1.000000)
         8       2. parameter for reaction( 2) : k8notP (k8notP_Reaction2) = 1000000.000000 (1000000.000000)
         9       3. parameter for reaction( 3) : k9 (k9_Reaction3) = 1000.000000 (1000.000000)
        10       4. parameter for reaction( 4) : k3 (k3_Reaction4) = 200.000000 (200.000000)
        11       5. parameter for reaction( 5) : k5notP (k5notP_Reaction5) = 0.000000 (0.000000)
        12       6. parameter for reaction( 6) : k1aa (k1aa_Reaction6) = 0.015000 (0.015000)
        13       7. parameter for reaction( 7) : k2 (k2_Reaction7) = 0.000000 (0.000000)
        14       8. parameter for reaction( 8) : k7 (k7_Reaction8) = 0.600000 (0.600000)
        15       9. parameter for reaction( 9) : k4 (k4_Reaction9) = 180.000000 (180.000000)
        16      10. parameter for reaction( 9) : k4prime (k4prime_Reaction9) = 0.018000 (0.018000)
        17
        18   --- Determine species role in reaction ------------------------
        19       1. species( 1) reaction( 1) : [0]
        20       2. species( 1) reaction( 2) : [0]
        21       3. species( 1) reaction( 3) : [0]
        22       4. species( 1) reaction( 4) : [0]
        23       5. species( 1) reaction( 5) : [0]
        24       6. species( 1) reaction( 6) : [0 1 0 0 1]
        25       7. species( 1) reaction( 7) : [1 0 0 1 0]
        26       8. species( 1) reaction( 8) : [1 0 0 1 0]
        27       9. species( 1) reaction( 9) : [0]
        28      10. species( 2) reaction( 1) : [1 0 0 1 0]
        29      11. species( 2) reaction( 2) : [0 1 0 0 1]
        30      12. species( 2) reaction( 3) : [1 0 0 1 0]
        31      13. species( 2) reaction( 4) : [0]
        32      ...
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga               ()            Brief Tutorial on Using DOTcvpSB             April 5, 2009   35 / 43
                                                                      SBML to DOTcvpSB      SBML to DOTcvpSB: Part 3



   SBML TO DOT CVP SB: PART 3

         1   --- Determine stoichiometry matrix       ----------------------------
         2       1. species(        EmptySet)         : 0 0 0 0 0 -1 1 1 0
         3       2. species(              C2)         : 1 -1 1 0 0 0 0 0 0
         4       3. species(              CP)         : 0 1 -1 -1 0 0 0 0 0
         5       4. species(               M)         : -1 0 0 0 -1 0 0 0 1
         6       5. species(              pM)         : 0 0 0 1 1 0 0 0 -1
         7       6. species(               Y)         : 0 0 0 -1 0 1 -1 0 0
         8       7. species(              YP)         : 1 0 0 0 0 0 0 -1 0
         9
        10   ---- Determine initial concentration --------------------------
        11       1. species(        EmptySet)   : 0.000000
        12       2. species(              C2)   : 0.000000
        13       3. species(              CP)   : 1.000000
        14       4. species(               M)   : 0.000000
        15       5. species(              pM)   : 0.300000
        16       6. species(               Y)   : 0.000000
        17       7. species(              YP)   : 0.000000
        18
        19   --- Get rate laws from reactions ------------------------------
        20       1. species(        EmptySet)   : 0
        21       2. species(              C2)   : + (k6_Reaction1*M) - (C2*k8notP_Reaction2) + (CP*k9_Reaction3)
        22       3. species(              CP)   : + (C2*k8notP_Reaction2) - (CP*k9_Reaction3) - (CP*k3_Reaction4*Y)
        23       4. species(               M)   : - (k6_Reaction1*M) - (k5notP_Reaction5*M) + (pM*(k4prime_Reaction9+k4_Reaction9*power(M,2)))
        24       5. species(              pM)   : + (CP*k3_Reaction4*Y) + (k5notP_Reaction5*M) - (pM*(k4prime_Reaction9+k4_Reaction9*power(M,2)))
        25       6. species(               Y)   : - (CP*k3_Reaction4*Y) + (k1aa_Reaction6) - (k2_Reaction7*Y)
        26       7. species(              YP)   : + (k6_Reaction1*M) - (k7_Reaction8*YP)
        27
        28   --- Get rate laws from rules ----------------------------------
        29       1. species(        EmptySet)   : 0
        30       2. species(              C2)   : 0
        31       3. species(              CP)   : 0
        32       4. species(               M)   : 0
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga         ()            Brief Tutorial on Using DOTcvpSB                              April 5, 2009   36 / 43
                                                                              SBML to DOTcvpSB      SBML to DOTcvpSB: Part 4



   SBML TO DOT CVP SB: PART 4

         1        5. species(                         pM)    : 0
         2        6. species(                          Y)    : 0
         3        7. species(                         YP)    : 0
         4
         5   --- Get species algebraic rules -------------------------------
         6       1. species(        EmptySet)   : 0
         7       2. species(              C2)   : 0
         8       3. species(              CP)   : 0
         9       4. species(               M)   : 0
        10       5. species(              pM)   : 0
        11       6. species(               Y)   : 0
        12       7. species(              YP)   : 0
        13
        14   --- Get species assignment rules               ------------------------------
        15       1. species(        EmptySet)               1.: 0
        16       2. species(              C2)               2.: 0
        17       3. species(              CP)               3.: 0
        18       4. species(               M)               4.: 0
        19       5. species(              pM)               5.: 0
        20       6. species(               Y)               6.: 0
        21       7. species(              YP)               7.: 0
        22
        23        The final order for replacement is: 1 7 5 3 2 6 4
        24
        25   --- Get information about events ------------------------------
        26
        27   --- Get information about rules -------------------------------
        28
        29   --- Get information about compartments ------------------------
        30       1. comp.(            cell)     : NaN
        31
        32        Please, select an action:
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                 ()            Brief Tutorial on Using DOTcvpSB             April 5, 2009   37 / 43
                                                                                SBML to DOTcvpSB      SBML to DOTcvpSB: Part 5



   SBML TO DOT CVP SB: PART 5

         1         [1] simulation of the model
         2         [2] single-optimization of the model
         3         1
         4
         5         Please, enter a final time (data.odes.tf) of the simulation: 100
         6
         7   ---   Save:   the original model from SBML              ...   OK
         8   ---   Save:   the information about the rules           ...   OK
         9   ---   Save:   the information about the species         ...   OK
        10   ---   Save:   parameters                                ...   OK
        11   ---   Save:   compartments                              ...   OK
        12   ---   Save:   events if they are define                 ...   OK
        13   ---   Save:   rules                                     ...   OK
        14   ---   Save:   differential equations                    ...   OK
        15   ---   Save:   the initial concentration                 ...   OK
        16
        17         Input file was saved correctly as dotcvp_Tyson1991_CellCycle_6var.m
        18         Please, check if the resulting input file is correct.
        19
        20   ans =
        21
        22                   typecode:         ’SBML_MODEL’
        23                     metaid:         ’_000001’
        24                      notes:         [1x3684 char]
        25                 annotation:         [1x2086 char]
        26                 SBML_level:         2
        27               SBML_version:         1
        28                       name:         ’Tyson1991_CellCycle_6var’
        29                         id:         ’Tyson1991CellModel_6’
        30         functionDefinition:         [1x0 struct]
        31             unitDefinition:         [1x0 struct]
        32                                     ...
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga             ()                  Brief Tutorial on Using DOTcvpSB             April 5, 2009   38 / 43
                                                                SBML to DOTcvpSB         SBML to DOTcvpSB: Simulation



   SBML TO DOT CVP SB: S IMULATION

         1   ________________________________________________________
         2
         3    DOTcvp: Dynamic Optimization Toolbox with CVP approach
         4     for handling continuous and mixed-integer DO problems
         5
         6     Main author: Tomas Hirmajer, thirmajer@iim.csic.es
         7                                  thirmajer@gmail.com
         8       Coauthors: Eva Balsa-Canto and Julio R. Banga
         9       Web pages: http://www.iim.csic.es/~dotcvp/
        10                  http://www.iim.csic.es/~dotcvpsb/
        11    Core version: DOTcvp_R2008_E2
        12   ________________________________________________________
        13
        14   ________________________________________________________
        15
        16    DOTcvp - a Module for Simulation
        17   ________________________________________________________
        18
        19   Generation of the file: cvm_rhs.m (ODE - MATLAB) .................. done!
        20   Simulation in progress, please wait ...............................
        21   Deleting of the temporary files ................................... done!
        22
        23   The detailed information is saved to the workspace structure with the name ’data’.




    Only the files, which are necessary for the simulation are generated and later deleted.

Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()               Brief Tutorial on Using DOTcvpSB                 April 5, 2009   39 / 43
                                                                                        SBML to DOTcvpSB      cdc2-cyclin Model



   G RAPHICAL O UTPUT: S TATE T RAJECTORIES
                                                                                           Tyson1991−CellCycle−6var
                                                                        1.4
                                                                                                                                             x
                                                                                                                                                 1
                                                                                                                                             x
                                                                                                                                                 2
                                                                        1.2                                                                  x
                                                                                                                                                 3
                                                                                                                                             x
                                                                                                                                                 4
                                                                                                                                             x
                                                                                                                                                 5
                                                                         1
                                                                                                                                             x
                                                                                                                                                 6
                                                                                                                                             x
                                                                                                                                                 7



                                                      State variables
                                                                        0.8



                                                                        0.6



                                                                        0.4



                                                                        0.2



                                                                         0
                                                                          0        20          40              60                 80   100
                                                                                                      Time


    F IGURE : The presented figure shows the dynamical behavior of the cdc2-cyclin model. The figure is generated by the simulation
    module.
Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                           ()            Brief Tutorial on Using DOTcvpSB                         April 5, 2009   40 / 43
                                                              SBML to DOTcvpSB      cdc2-cyclin Model



   T HANK YOU FOR YOUR ATTENTION




                                                           Thank you for your Attention




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga   ()          Brief Tutorial on Using DOTcvpSB      April 5, 2009   41 / 43
                                                                            Toolbox References



   T OOLBOX R EFERENCES

                 IVP S OLVER :
                    SUNDIALS A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward, Sundials: Suite of nonlinear and
                             differential/algebraic equation solvers, ACM Transactions on Mathematical Software, 31 (2005), pp. 363-396.
               NLP S OLVERS :
                          IPOPT A. Wächter and L. T. Biegler, On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear
                                programming, Mathematical Programming, 106 (2006), pp. 25-57.
                     FMINCON T. Coleman, M. A. Branch, and A. Grace, Optimization toolbox for use with matlab user’s guide ver. 2, 1998.
                               DE R. Storn and K. Price, Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces, Journal of
                                  Global Optimization, 11 (1997), pp. 341-359.
                            SRES T. P. Runarsson and X. Yao, Stochastic ranking for constrained evolutionary optimization, IEEE Transactions Evolutionary Computation, 4
                                 (2000), pp. 284-294.
           MINLP S OLVERS :
                          MISQP O. Exler and K. Schittkowski, A trust region sqp algorithm for mixed-integer nonlinear programming, Optimization Letters, 1 (2007), pp.
                                269-280.
                          ACO MI M. Schlüter, J. A. Egea, and J. R. Banga. Extended ant colony optimization for non-convex mixed integer nonlinear programming.
                                 Computers & Operations Research, (2008), doi:10.1016/j.cor.2008.08.015.
                            MITS O. Exler, L. T. Antelo, J. A. Egea, A. A. Alonso, and J. R. Banga, A tabu search-based algorithm for mixed-integer nonlinear problems and
                                 its application to integrated process and control system design, Computers and Chemical Engineering, 32 (2008), pp. 1877-1891.




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga            ()               Brief Tutorial on Using DOTcvpSB                                               April 5, 2009   42 / 43
                                                              Problem References and Acknowledgments



   P ROBLEM R EFERENCES AND ACKNOWLEDGMENTS



            S INGLE O PTIMIZATION :
    DRUG DISPLACEMENT PROBLEM                 E. Balsa-Canto, V. S. Vassiliadis, and J. R. Banga. Dynamic optimization of single- and multi-stage systems using a hybrid
                                              stochastic-deterministic method. Industrial and Engineering Chemistry Research, 44(5):1514-1523, 2005
      CALCIUM OSCILLATOR MODEL                U. Kummer, L. F. Olsen, C. J. Dixon, A. K. Green, E. Bornberg-Bauer, and G. Baier. Switching from Simple to Complex Oscillations
                                              in Calcium Signaling. Biophysical Journal, 79(3):1188-1195, 2000.
   CALCIUM OSCILLATOR PROBLEM                 Lebiedz, D., Sager, S., Bock, H., and Lebiedz, P. (2005). Annihilation of limit-cycle oscillations by identification of critical perturbing
                                              stimuli via mixed-integer optimal control. Phys. Rev. Lett., 95, 108303.
             SBML TO DOT CVP SB:
                 CDC 2- CYCLIN MODEL          J. J. Tyson, Modeling the cell division cycle: cdc2 and cyclin interactions, Proc. Nati. Acad. Sci. USA, 88 (1991), pp. 7328-7332.
                       SBMLT OOLBOX S. M. Keating, B. J. Bornstein, A. Finney, and M. Hucka, SBMLToolbox: an SBML toolbox for MATLAB users, Bioinformatics, 22
                                    (2006), pp. 1275-1277.
               ACKNOWLEDGMENTS : The first author kindly acknowledges the support under the JAE-doc programme from the CSIC, Spain.




Tomáš Hirmajer, Eva Balsa-Canto, and Julio R. Banga                ()                  Brief Tutorial on Using DOTcvpSB                                                April 5, 2009   43 / 43

								
To top