How to set up steady state constraints

Preface

Often, the dynamics of the model should be in a steady state at the beginning of a simulation. This reflects the biological situation of cells being in a resting (i.e. unstimulated) state. D2D offers three approaches to deal with steady states; each of which have their own advantages and drawbacks:

Analytical steady states: Analytical steady states via parameter transformations in the CONDITIONS field

The steady state of an ODE system \dot{x} = f(x) can be calculated analytically by solving 0 = f(x). The solution can be implemented by specyfing the respective constrains the CONDITIONS section of the model or data .def file. When the analytical steady state is known, its usage is superior to pre-equiblibrating the system as it reduces the simulation workload. However, solving the ODE system by hand is only feasible for very simple systems.

The python tool ODESS (Rosenblatt et al., 2016) is available to analytically determine steady states of ODE systems. To use it from d2d, the model must be exported as a stochiometric matrix in csv format, which can be achieved by running arExportModelToDmod after compilation.

Pre-equilibration based steady states: Steady state imposition by equilibrating the system to steady state.

Advantages: No additional term in optimization. Reduction in the number of parameters.

Disadvantages: Pre-simulation required. Integration failure when steady state does not exist.

This method is described here

Constraint based steady states: Steady state imposition by means of penalized optimization.

Advantages: Optimized initial condition parameters correspond to a (near) steady state and it is also possible to consider systems 'near' steady state.

Disadvantages: Constraint is included in profiling and optimization, which can trade off penalty and fidelity to the data. Have to critically evaluate what the parameter profiles really mean when introducing such a constraint. Since constraints are ratio-metric, states near zero are penalized much more strongly in absolute terms.

This method is described on this page

Using constraint based steady states

A steady state constraint can be assigned by invoking:

ar.model(jm).condition(jc).qSteadyState(:) = 1;

for each model with index jm and experimental condition with index jc. It is important to figure out for which experimental condition jc the steady state constraints should be activated. Here the function arShowDataConditionStructure can help. It is also important to make sure that the specified condition can really lead to a steady state.

The violation of the constraints are then displayed in the output of the arChi2 function, and are used in the optimization of the model.

The constraints are implemented as soft constraints. The strength of the constraint to be fulfilled is controlled by the field:

ar.model(jm).condition(jc).stdSteadyState(:) = std_val;

It indicates the deviation from zero of the first derivative of the dynamics with respect to time dx/dt at the beginning of the simulation of the respective experimental condition. The smaller stdSteadyState, here given by std_val, the stronger the constraint. In principle each constraint component can have an individual stdSteadyState value, however, this is probably not necessary in most cases.

The remaining question is of course, what a reasonable value for std_val is. Here, we can give a hint: Fulfilling the constraint is an artificially imposed, i.e. not measured, requirement to the model. It is reasonable to weight it equally in the optimization than the entirety of experimental data. Also, it is important to consider the typical timescale that is of interest. For instance, the cells might be allowed to change over time very slowly. Here is an example script that selects the stdSteadyState value accordingly:

% turn on steady state constraint
ar.model(1).condition(1).qSteadyState(:) = 1; arChi2;

% allow approx. 1e-1 relative deviation from steady state with respect to 200 minutes simulation of dynamics
ar.model(1).condition(1).stdSteadyState(:) = 1e-1 / 200;

% equally weight constraints versus experimental data
ar.model(1).condition(1).stdSteadyState = ar.model(1).condition(1).stdSteadyState * (ar.nconstr/ar.ndata);