Data2Dynamics / d2d

a modeling environment tailored to parameter estimation in dynamical systems
https://github.com/Data2Dynamics/d2d
57 stars 29 forks source link

How to include conservation relation constraints? #103

Closed SchSascha closed 6 years ago

SchSascha commented 6 years ago

Hi,

I want to include conservation relations and wonder how I can make sure the relation holds, such that the total species concentration is always bigger than its variants. Example: say, we know X1 and X2 together form X_total by X_total = X1(t) + X2(t) and follow X1(t)/dt = b*X2(t) - a*X1(t) X2(t)/dt = a*X1(t) - b*X2(t) I want to reduce it to X2(t)/dt = a*(X_total - X2(t)) - b*X2(t) with X_total being a free parameter and the subtraction (X_total - X2(t)) resembling X1(t) quantities. However, this subtraction can result in negative values upon fitting data, which is obviously unwanted. Is there a way to include parameter constraints like X_total >= X2(t) for all time points in the setup of the model file?

Best Sascha

JoepVanlier commented 6 years ago

There might be ways to introduce it as a custom residual or in the form of a data and observable pair, but I would not recommend this route, as it may cause other problems down the road.

If it's a conserved moiety in the system, wouldn't it already always be fulfilled as long as you make sure that your initials are chosen correctly?

As for the initials, you could express the initials for the pools as a fraction of the total in the CONDITIONS field and constrain that parameter between zero and one.

Just out of curiosity, is there a specific reason you want to eliminate the conserved moieties? There is already a function that does this for the model internally (see arReduce). It will still have separate initials for all the species involved in the conserved moiety, but internally it handles the substitutions.

SchSascha commented 6 years ago

Hi Joep,

thanks for your fast reply!

The reason I want to include it myself is that I have several datasets, including different perturbations that are associated to the same model. However, at time point zero, some datasets have already seen inhibitors, hence require different initial concentrations. As I have semi-quantitative data, I do not know the total concentrations, but want them to be calibrated when fitting the model. As I want to test neglecting protein synthesis/degradation I want to make sure that E_total is the same over all datasets. I thought one convenient way to achieve this is via conservation pools/relations, as then the E_totals are free parameters, but get calibrated over all datasets (hence the same E_totals over all relevant species). I would believe using arReduce (thanks for the hint!) would not do the job, would it? Expressing the initials as fractions of the total sounds an interesting approach. My first thought then is how can I make sure that the sum of all initial fractions that are associated to a specific E_total equal one?

In the end my question about conservation pools really is just a try to tackle the question about how to configure that the sum of (relative) concentrations of a subset of species in the model stays the same (aka gets calibrated to the same arbitrary concentration) across several datasets?

JoepVanlier commented 6 years ago

Hey Sascha,

I would split out the initials into fractions. That way you avoid having to explicitly impose a total concentration relation.

With one it is simple, then the CONDITIONS field would look something like this for your example: init_stateA "f_A * init_total" init_stateB "(1-f_A) * init_total"

It is not necessary to remove the total pool in this case either. Only if you want to for some speific reason. Note that if f_A is dataset specific, you might want to call it f_A_dataset1 or something. If you know it will be different for every dataset that you load, then you might even write it as f_A_filename in the model file and not worry about it again. Filename then gets replaced by the name of the def file.

If you need to divide the total over more states, you just repeat the trick. Single split, the total would be: fa*t + (1-fa)*t

We can split one of these states up again to get: fa*t + (1-fa)*fb*t + (1-fa)*(1-fb)*t

This sums back to the same total and as long as you choose fa and fb between 0 and 1 (via normal independent parameter bounds) this should keep all states positive. The individual initials would then be fa*total (1-fa)*fb*total (1-fa)*(1-fb)*total

Best, Joep.

SchSascha commented 6 years ago

Hey Joep,

this sounds really nice, haven't thought about it this way. I have up to four states in one pool, but the math is easy here.

Will give this a go, thank you!

Best Sascha

JoepVanlier commented 6 years ago

Glad to help. Feel free to reopen the issue if you run into trouble with it :)