steady state fitting - set up

[carlotta]

Dear all,

I am trying to do a steady-state fitting of my data but I am not sure I am doing it correctly since the model is not fitting (see image attached). 

I am mixing two species (A and B) in different amounts which complexate in the species AB.

AB and B are then in equilibrium with two other species: C and D. This mechanism I think can be written with the following two equilibria reactions:

 A + B =  AB

AB + B = C + D

 Since I am used to time-dependent data with Copasi, I am not sure how should I do with this set of time-independent data, and I would have some questions:

-  I would go for Model  -> Biochemical and write the two equations. But then what should I set as initial concentration for the species? The concentration at the equilibrium (which is the one that I measured) or the concentration of A and B that I used to prepare the reactions?

 - In the Task -> Parameter estimation, I upload my experimental data. Which are my independent values?

- Is there something else I should consider since is a steady-state fitting?

Thank you a lot for the help

Best regards

Carlotta

[Frank Bergmann]

Hello Carlotta, 

To your questions: 

 - i'd use the concentration of A and B. 

- as for the independent data. Are you sure you want to choose the InitialParticleNumber of the species? Maybe the initial concentrations would be a better choice. But i cant say for sure without looking at the data. 

- I dont think there is anything else specifically to consider that would differentiate time course parameter estimation from steady state parameter estimation in COPASI. The one thing to keep in mind for it generally is that it will only work of the model can reach a steady state. But in your case that is no problem, as you will reach an equilibrium steady state.  

I'd be happy to take a look if you were to send me the model. 

best

Frank

--
You received this message because you are subscribed to the Google Groups "COPASI User Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to copasi-user-fo...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/e519130c-7e34-4c4d-ae05-6f9e4130b2fbn%40googlegroups.com.

[carlotta]

Dear Frank,

Thanks a lot for the answer.

- as initial concentration for A and B I set the initial amount, not the one at equilibrium

- yes, sure my mistake, it is better to choose the initial concentrations

Please find attached the model and the .txt file with the data.

Best

Carlotta

[Frank Bergmann]

Hello Carlotta, 

you might not have gotten the mail i sent to you privately before. 

Your parameter estimation setup looks fine to me, we have 

- C, D, AB mapped as dependent species, indicating that these were the measured values at the end of the experiment

- we have A and B marked as independent species, indicating that these are the concentration provided for the experiment (do we know what the initial concentrations for C and D are at the start of the experiment)

The only thing i would change is the boundaries for the parameters to vary. Running parameter estimation reaches those very quickly so it would make sense to increase them. 

I think the reason for the bad fit, is that the model structure does not support the data. 

If we have the system of the two reversible mass action reactions: A + B = AB and AB + B = C + D. Then I would expect that at steady state when we have C and D present, that there would be some AB as well. 

But maybe someone else has a better idea. 

best

Frank

[carlotta]

Dear Frank,

I have figured out how to fit my data with copasi and now the model is working!

I would have still a couple of questions.

 

Here a short recap of the project.

 I am mixing two species (B and T) in different amounts five different times. Each of these gives a condition corresponding to a row in my .txt file.

In the first two, the amount of B is bigger than T, in the third B and T are the same, while in the last two T is the species in excess.

The system can be described with two equilibrium reactions:

B + 2 * T = BT2

B + BT2 = AT2 + C

 

I have added two columns to my old .txt file: one with the initial values of B (B_0) and the other with the initial values of T (T_0),  the amounts that I mixed initially. With this addition, the program is now fitting nicely.

 

So in my .txt file there are 7 columns, all expressed in concentration (mol/L): B_0, T_0, B, T, BT2, AT2, C. And 5 rows.

Since B_0 and T_0 are not described by the reactions, I added them as species with a simulation type “assignment”. The mathematical expression I used to describe B_0 is the initial concentration of B ({[B]_0}) while T_0 is the initial concentration of T ({[T]_0}).

 

To calculate I used the Parameter Estimation task, where I assigned the 4 constants k to the two equilibria (forward and backward reactions).

In the Experimental Data, I loaded the .txt file, and defined B_0 and T_0 as independent values, while all the others as dependent. The experiment type is steady-state.

 

When I run the program the data are fitted very well, but there are some problems with the values of the rate constants.

-          Under Parameter Estimation --> Results --> Parameters, the coefficients of variation are very high (1e+8 %). How important is it that these have reasonable values? What could I change in the model to make them more reliable?

-          I am rather interested in the ratio between the rate constants since it is a thermodynamic state. But every time I run the program, it gives me different values for the k, and a different ratio. Maybe this is related to the previous problem.

 

Attached is the program the .txt file.

Thank you a lot for all your help

Best regards

Carlotta

[sven....@gmail.com]

Hello Carlotta, 

this is a question of structural parameter identifiability (or rather unidentifiability in this case). Basically the very high coefficients of variation indicate that it is not possible to determine the value of all parameters at the same time from the available data. 

This means that different combinations of parameter values would lead to exactly the same fit. 

There is a more mathematical definition of this, but in this case we can understand directly why this is the case: All unknown parameters are rate constants (that describe the speed of a reaction) but the data, since it is steady state data, has no information at all about any time scale. Therefore the time scale of the model is completely undetermined. This is also the reason why you can get completely different parameter values if you repeat the estimation (but the fit looks good each time). 

You already have the solution to that: The ratio of forward and backwards rate constants should be identifiable, since it does not depend on the time scale. There is two ways to implement that: You could rewrite the rate law for the reactions, so that it explicitly uses the equilibrium constant, and if you then rerun the parameter estimation you should get a more reasonable coefficient of variation for those. 

The easier method is to just fix either the forward or the backwards rate constant (of both reactions) to an arbitrary value (e.g. 1.0). By doing this you basically provide an (arbitrary) time scale to the model, and the remaining 2 parameters should then be identifiable. The fact that the timescale of the model is now fixed arbitrarily should not be a problem as long as you only use the model to explain steady state experimental data. 

I have tried the easier method (using your model) and it works for me. 

However, I still get a coefficient of variation of more than 100% for the first reaction's parameter (the second reaction is fine with 8%). This is yet another effect that is called "practical nonidentifiability". That means a parameter can be identified in principle, but with such a low accuracy that it may be useless.  I think I know the reason for that: In the data for most experimental conditions the value for T is exactly 0.0. That means in order to get a good fit COPASI will choose an extremely large value for the forward rate constant of the first reaction to make sure that the equilibrium is very much on the BT2 side of the reaction and T is very small. This is correct and understandable. However, the exact value for this very large rate constant is not reliably identifiable. Any very large value can lead to a good fit, and you can change the value of the parameter quite substantially without the fit getting much worse. 

There is not much you can do about this with the current data set, basically the qualitative result is "the equilibrium of the  first reaction is very much towards BT2"

Hope that helps, 

Sven

[carlotta]

Dear Sven,

Thanks a lot for the answer and the analysis, it is very useful!

I just have some questions about your comments:

- You said <<You could rewrite the rate law for the reactions so that it explicitly uses the equilibrium constant, and if you then rerun the parameter estimation you should get a more reasonable coefficient of variation for those>>

How can I define the equilibrium constant? I thought to create a global quantity K2 and chose as type "assignment" giving then the mathematical expression k1(2) / k2(2), and in parameter estimation, I would not add k1(1) and k2(1) as parameters. Is this the right way? Where can I then find the value of K2 after running the fitting?

- How can I fix a rate constant to an arbitrary value?

- If then the equilibrium of the first reaction is very much towards BT2, would it make more sense to define it as an irreversible reaction rather than an equilibrium with a very high equilibrium constant?

Thank you

Best regards

Carlotta

[sven....@gmail.com]

Dear Carlotta, 

carlotta schrieb am Donnerstag, 23. Februar 2023 um 13:56:21 UTC+1:

Dear Sven,

Thanks a lot for the answer and the analysis, it is very useful!

I just have some questions about your comments:

- You said <<You could rewrite the rate law for the reactions so that it explicitly uses the equilibrium constant, and if you then rerun the parameter estimation you should get a more reasonable coefficient of variation for those>>

How can I define the equilibrium constant? I thought to create a global quantity K2 and chose as type "assignment" giving then the mathematical expression k1(2) / k2(2), and in parameter estimation, I would not add k1(1) and k2(1) as parameters. Is this the right way? Where can I then find the value of K2 after running the fitting?

This would be a possible way to do it (if I understand you correctly): 

-Create a global quantity Keq (that you add to the parameter estimation procedure)

-Create a global quantities k1 and k2, and define an assignment for k2 (k2 := k1 / Keq)

-Set an arbitrary value for k1 in the model and do not add k1 and k2 to the parameter estimation. 

A different approach would be to create a rate law that includes Keq.  In your model for reaction 1 this could be something like

v = k ( S1*S2*S2 - 1/Keq*P)

where when you use it, S1 will be B, S2 will be T, and P will be BT2. If you then do a parameter estimation for k and Keq, the result will probably be that k is unidentifiable, but Keq may be identifiable. However, since this approach is more work than your suggestion above, I would only do this in a larger model. 

 

- How can I fix a rate constant to an arbitrary value?

 by fixing the rate constant I just meant not including it in the fitting procedure. Then the value will just stay as it is defined in the model (e.g. 1.0). 

 

- If then the equilibrium of the first reaction is very much towards BT2, would it make more sense to define it as an irreversible reaction rather than an equilibrium with a very high equilibrium constant?

That depends on the biological (or biochemical) context. With the current data and if you think that it is correct that T is consumed completely whenever there is enough B, then reaction 1 can be modelled as irreversible. If, on the other hand, you expect to get more data later where T is not exactly 0.0 but rather a small value, or if you plan to use the model for other circumstances, then it might make sense to keep it reversible.

Sven

 

[carlotta]

Thanks a lot Sven, it was very useful!

Best

Carlotta

[carlotta]

Dear Sven, 

I have another set of steady-state data that I would like to fit with equations similar to the last time.

I am mixing two species (C and T) in different amounts five different times. Each of these gives a condition corresponding to a row in my .txt file.

The system can be described with three equilibrium reactions:

C + 2 * T = CT2

C + CT2 = BT2 + D

BT2 + D = AT2 + E

 In my .txt file I have two columns with the initial values of C(C_0) and T(T_0), which are the amounts that I mixed in the beginning. In addition, I can measure the amount of AT2, BT2, CT2 and T at the equilibrium, and for each of these species, I added a column in the .txt file. I don’t know the values of C, D and E, so they will not have a column in the .txt. 

       As previously suggested by you, I fix the backreactions to 1, and in the parameter estimation task, I only calculate the three forward reactions. In the Experimental Data, I loaded the .txt file, defining C_0 and T_0 as independent values, while all the others as dependent. The experiment type is steady-state.

       But I have some doubts regarding the fittings:

       - When I run the program the data are not fitted well and I am not sure why, since for the other case it worked. (In the previous case, for all the species present in the reactions a column in the .txt file was assigned. But I also tried to simulate the system excluding some variables and it works nicely anyway.)

        - Since I don’t know the value of all the species present, is there a way to simulate their concentrations? I know that in the case of a kinetic simulation, after calculating the rate constants, the time course can be simulated and all the species can be plotted vs time. Is there a similar simulation for the steady state? I tried to simulate it with the task “steady state”, but as outcome, I only get one value for each species.

-          - I realized that we are forcing a program thought for fitting kinetic to fit the equilibrium data. So I am not sure Copasi is the right software to use for this project. Do you maybe know something more suitable? Or do you think copasi is indeed suitable also for these kind of tasks?

Attached is the program and the .txt file.

Thank you a lot for all your help

Best regards

Carlotta

[sven....@gmail.com]

Dear Carlotta,

I think the setup of the parameter estimation is correct (as far as Copasi-technicalities are concerned).

carlotta schrieb am Donnerstag, 6. April 2023 um 01:17:16 UTC+9:

I have another set of steady-state data that I would like to fit with equations similar to the last time.

I am mixing two species (C and T) in different amounts five different times. Each of these gives a condition corresponding to a row in my .txt file.

The system can be described with three equilibrium reactions:

C + 2 * T = CT2

C + CT2 = BT2 + D

BT2 + D = AT2 + E

 In my .txt file I have two columns with the initial values of C(C_0) and T(T_0), which are the amounts that I mixed in the beginning. In addition, I can measure the amount of AT2, BT2, CT2 and T at the equilibrium, and for each of these species, I added a column in the .txt file. I don’t know the values of C, D and E, so they will not have a column in the .txt. 

       As previously suggested by you, I fix the backreactions to 1, and in the parameter estimation task, I only calculate the three forward reactions. In the Experimental Data, I loaded the .txt file, defining C_0 and T_0 as independent values, while all the others as dependent. The experiment type is steady-state.

This should be fine.

 

       But I have some doubts regarding the fittings:

       - When I run the program the data are not fitted well and I am not sure why, since for the other case it worked. (In the previous case, for all the species present in the reactions a column in the .txt file was assigned. But I also tried to simulate the system excluding some variables and it works nicely anyway.)

It is definitely possible to do a parameter estimation when only some of the variables of the model are measured.

        - Since I don’t know the value of all the species present, is there a way to simulate their concentrations? I know that in the case of a kinetic simulation, after calculating the rate constants, the time course can be simulated and all the species can be plotted vs time. Is there a similar simulation for the steady state? I tried to simulate it with the task “steady state”, but as outcome, I only get one value for each species.

If you run the Steady State task (and if it succeeds) you will get a table with the steady state concentrations of all the variables.  If you want to see the steady states for different starting values of C and T, you have to change these values in the model (or do a parameter scan, see below). 

I do not quite understand what you mean by "simulation". You could say the steady state task is the simulation of the steady state behaviour of the system. Of course you can also do a time course simulation of the model, but the dynamic behaviour is not really meaningful since you have only fitted the model to steady state data. 

-          - I realized that we are forcing a program thought for fitting kinetic to fit the equilibrium data. So I am not sure Copasi is the right software to use for this project. Do you maybe know something more suitable? Or do you think copasi is indeed suitable also for these kind of tasks?

This is perfectly fine. COPASI is definitely also designed to work with steady state data. 

As far as I can see the fact that you cannot get a good fit is actually the correct result, meaning that the model as it is now cannot fit the data. This is actually a valid result of a parameter estimation, and it tells you  that either the model is missing something, or that something is wrong with the data. I have no idea which of the two is the case. 

 

What I can see from the model is the following: 

 

BT2 is measured, and it has its maximal value  in the "middle" experiment, where C0 and T0 are both 0.125. The model captures this behaviour nicely.

BT2 will alway be equal to D (since they are produced and consumed always in parallel).

Since AT2 is produced exclusively from BT2 and D (in the model), it must also have its maximum value under the same conditions, i.e. when C0 and T0 are at 0.125. This is obviously contradicted by the data. What than means is that the model is absolutely not capable of reproducing the measurements for BT2 and AT2 at the same time, no matter the parameter values. 

Hope this helps, 

Sven

[carlotta]

Dear Sven,

Thanks a lot for all your help! I was able to fit my steady-state data with Copasi as you suggested, but without your feedback, it would have been very difficult.

Thanks again

Best

Carlotta