Problems with accessing dv/dS (Jacobian) #910
Comments
|
If it's dv/ds you want (where v is a reaction rate, not a rate of change)
then the Jacobian won't help. The Jacobian is defined as (dsi/dt)/dsj. It
is true that the jacobian is computed from sums of dv/ds but to get the
derivatives, dvds you need to call the get unscaled elasticities, ie
r.getuEE ("v", "s")
The numerical error you got is a problem, can you send us the model to look
at? I think what it means is that one of the species is close to zero and
its finding it difficult to compute dv/ds
Hebrert
…On Mon, Nov 8, 2021 at 7:05 AM Matthias König ***@***.***> wrote:
Hi all,
I want to access the changes in reaction rates depending on species.
As I remember this is the information available via the Jacobian, i.e.
dv_i/dS_k is the J_ik.
[image: image]
<https://user-images.githubusercontent.com/900538/140765029-52643266-decf-44be-adbb-2600762afb81.png>
So I tried to access the reduced Jabobian, but get a matrix which has
species ids for rows and columns (and not rate ids for rows and species ids
for columns)?
apap_ext, apap, napqi, CYP2E1
apap_ext [[ -0.0284312, 0.0927807, 0, 0],
apap [ 0.0284312, -0.099167, 0, -3.40902e-05],
napqi [ 0, 0.00638628, -0.000799083, 3.40902e-05],
CYP2E1 [ 0, 0, 0, 0]]
Is this correct? I would expect reaction ids in the rows here? What do the
entries mean in this square matrix? And how would I access the changes of
reaction rates with species?
When trying to access the full Jacobian, I get
python: /__w/1/s/source/rrRoadRunner.cpp:3765: double rr::RoadRunner::getUnscaledSpeciesElasticity(int, int): Assertion `std::abs(originalConc - tmp) < 1e-13' failed.
Process finished with exit code 134 (interrupted by signal 6: SIGABRT)
so that I cannot access the values.
Code and model attached:
apap_core.zip
<https://github.com/sys-bio/roadrunner/files/7497905/apap_core.zip>
import roadrunner
from pkdb_models.models.acetaminophen_pm import MODEL_BASE_PATH
from pkdb_models.models.acetaminophen_pm.models import apap_liver_core
# -------------------------------------------------------------------------------------
# load model
# -------------------------------------------------------------------------------------
# model_path = MODEL_BASE_PATH / f"{apap_liver_core.mid}_{apap_liver_core.version}.xml"
# change this to the following
model_path = "apap_liver_core_9.xml"
print(model_path)
# roadrunner instance for local point
r = roadrunner.RoadRunner(str(model_path))
# -------------------------------------------------------------------------------------
# selections
# -------------------------------------------------------------------------------------
# in a first step it must be defined what should be simulated in the model.
# The following sets this selection
r.timeCourseSelections = [
"time",
# local concentrations
"[apap_ext]",
"[apap]",
"[napqi]",
# sinks/sources
"d_apap_ext__dt",
# tangents
"d_vapap_ext__d_apap_ext",
# other readouts
"necrosis",
"CYP2E1",
# rate of change
"apap_ext'", # [mmole/min]
"apap'", # [mmole/min]
"napqi'", # [mmole/min]
]
# -------------------------------------------------------------------------------------
# set integrator settings
# -------------------------------------------------------------------------------------
# Some integrator settings are required to handle the very small volumes.
# set tolerances for very small FEM volumes
integrator: roadrunner.Integrator = r.integrator
integrator.setValue("absolute_tolerance", 1e-14)
integrator.setValue("relative_tolerance", 1e-14)
# -------------------------------------------------------------------------------------
# set values
# -------------------------------------------------------------------------------------
# volumes based on FEM point
# Vext is fluid phase of the point, Vcell is the fat + cell phase
# e.g. if you have
vol_point = 0.1 # [liter]
f_fluid = 0.2 # [dimensionless]
f_fat = 0.3 # [dimensionless]
f_tissue = 1 - f_fluid - f_fat
# setting volumes
r.setValue('Vext', vol_point * f_fluid) # [liter]
r.setValue('Vli', vol_point * (f_fat + f_tissue)) # [liter]
# protein amount based on position
# (varied with position of grid point between 0.0 periportal and 1.0 pericentral)
r.setValue('CYP2E1', 1.0) # [dimensionless] pericentral point
# concentrations of apap in fluid phase (local concentration of you fluid phase)
r.setValue('[apap_ext]', 1.0) # [mM] = [mmole/liter]
# -------------------------------------------------------------------------------------
# simulate timestep
# -------------------------------------------------------------------------------------
time = 0.0 # [min]
delta_time = 0.1 # [min]
# simulate a step
s = r.simulate(start=time, end=time+delta_time, steps=1)
print(s)
# -------------------------------------------------------------------------------------
# access results
# -------------------------------------------------------------------------------------
# source and sink terms
print("d_apap_ext__dt", s["d_apap_ext__dt"], "[mmole/min]")
# tangents dv/dc = [mmole/min]/[mmole/liter]
print("d_vapap_ext__d_apap_ext", s["d_vapap_ext__d_apap_ext"], "[l/mmole]")
# necrosis
print("necrosis", s["necrosis"], "[dimensionless]")
# internal concentrations
print("[napqi]", s["[napqi]"], "[mM]")
print("[apap]", s["[apap]"], "[mM]")
Jred = r.getReducedJacobian()
print(Jred)
J = r.getFullJacobian()
print(J)
—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
<#910>, or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAIBSDVZNOEGMAWRAX4WCFDUK7RDRANCNFSM5HS4MSZA>
.
Triage notifications on the go with GitHub Mobile for iOS
<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675>
or Android
<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.
--
Herbert Sauro, Professor
University of Washington, Bioengineering
206-685-2119, www.sys-bio.org
***@***.***
Books: http://books.analogmachine.org/
|
|
@hsauro Thanks. The This should be working from the description in The model was attached in the original comment. |
|
These selections should work, try lower case uee and ee. If that doesn't
work the current work around is to use oneStep.
|
|
I did not reproduce your errors with your models on my machine. My results were: I also discovered that to get the unscaled elasticity coefficient, the syntax is:
However, if I add this to your time course selections list in the program you provided, it fails, giving the following error: "The second argument to a metabolic control coefficient selection, apap_ext, must be either a global parameter, boundary species or conserved sum" Given that 'getuEE()' works with those arguments, this seems like a bug in roadrunner to me, perhaps having to do with it checking whether it's a metabolic control coefficient instead of an unscaled elasticity coefficient? At any rate, I'll look into it. |
|
Yes it should be able to compute the elasticity under any condition.
|
|
@luciansmith Thanks for looking into this. Did you change the volumes as in the scripts before calculating the Jacobian? If I just run the model with the default volumes I do not get the errors in the Jacobian. In addition I set stricter tolerances for smaller volumes ( It is strange to get an error |
|
So! A few things:
The fixes are incorporated into #911 Here's the updated python script I used, and the results I get from it: And the results: |
|
@luciansmith
I work with the libroadrunner and libroadrunner-experimental packages or download arzure artifacts for the develop branch (roadrunner-ManyLinux2014-py38 or py39). I never build from scratch. So the only way the debug mode is selected is if it is selected in your pip releases for linux (or the azure pipeline). Why are the reduced Jacobian and Jacobian different? Or is the reduced Jacobian something different? |
|
The top rows of the reduce jacobian should correspond to the same top rows
in the full jacobian, there is something wrong here at least in this model.
I am assuming you computed both of these at exactly the same state (doesn't
have to be the steady state)?
This will require a detailed look, Lucian if you have some time we could do
a zoom debugging session if that would help. just need to set up a time.
|
|
I noticed that the reduced and full are the same dimension so there
hasn't been any reduction (which is ok). In that case the two Jacobians
should be exactly the same.
I'll try out some models myself.
|
|
@hsauro @luciansmith So it seems the Azure wheels are build with debug mode on. This should probably be changed. |
|
Aha! OK, looks like our azure build didn't set the build type to 'release' for the manylinux builds. Thank you! However, now I'm wondering if this was fortuitous, because it indicates that the sanity checks have actually failed. I'll investigate further. @hsauro , if I need help somewhere I'll ping you and set up a zoom call. Thanks! |
|
Also, @hsauro : as you check things, try changing some of the values with 'setValue' before calculating the Jacobians. My guess is that one of our routines doesn't take that into account, somehow. |
|
Oh, and @matthiaskoenig : are you OK with your model going into our code as a test? I guess it's already available publicly here on this issue, but that's a bit more obscure than source code. |
|
ok I'll try out something now.
H
…On Tue, Nov 9, 2021 at 3:27 PM Lucian Smith ***@***.***> wrote:
Also, @hsauro <https://github.com/hsauro> : as you check things, try
changing some of the values with 'setValue' before calculating the
Jacobians. My guess is that one of our routines doesn't take that into
account, somehow.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
<#910 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AAIBSDVRFLPSRIN5MH377Z3ULGU4XANCNFSM5HS4MSZA>
.
Triage notifications on the go with GitHub Mobile for iOS
<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675>
or Android
<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.
--
Herbert Sauro, Professor
University of Washington, Bioengineering
206-685-2119, www.sys-bio.org
***@***.***
Books: http://books.analogmachine.org/
|
|
First of all a basic sanity check I generated this random model that has a
stead state and compute the two jacobian, as you can see they are the same.
I did a bunch of setValues and no difference. I'll now look at Matthias'
model
r = te.loada ('''
var S0, S1, S2, S5, S6, S8, S9
ext S3, S4, S7;
J0: S0 -> S6; E0*(k0*S0);
J1: S9 + S2 -> S6; E1*(k1*S9*S2);
J2: S2 -> S7; E2*(k2*S2);
J3: S6 -> S7 + S5; E3*(k3*S6);
J4: S9 + S8 -> S6; E4*(k4*S9*S8);
J5: S3 -> S0; E5*(k5*S3);
J6: S6 -> S1; E6*(k6*S6);
J7: S6 -> S4 + S0; E7*(k7*S6);
J8: S5 -> S9; E8*(k8*S5);
J9: S2 -> S9 + S5; E9*(k9*S2);
J10: S1 -> S8; E10*(k10*S1);
J11: S3 + S1 -> S8; E11*(k11*S3*S1);
J12: S9 -> S2; E12*(k12*S9);
J13: S9 -> S4 + S7; E13*(k13*S9);
J14: S2 -> S9 + S0; E14*(k14*S2);
r.setValue ('S2', 1)
r.setValue ('S8', 1)
r.setValue ('S9', 1)
r.setValue ('S1', 1)
r.setValue ('S5', 1)
S0, S1, S2, S5, S6, S8,
S9
S0 [[ -0.41782, 0, 0.212978, 0, 0.8178, 0,
0],
S1 [ 0, -0.772176, 0, 0, 0.0419199, 0,
0],
S2 [ 0, 0, -1.45948, 0, 0, 0,
0.7046],
S5 [ 0, 0, 0.438294, -0.559132, 0.404227, 0,
0],
S6 [ 0.41782, 0, 0.330554, 0, -1.26395, 0.616207,
0.346726],
S8 [ 0, 0.772176, 0, 0, 0, -0.616207,
-0.140416],
S9 [ 0, 0, 0.320719, 0.559132, 0, -0.616207,
-2.03405]]
r.getReducedJacobian()
Out[521]:
S0, S1, S2, S5, S6, S8,
S9
S0 [[ -0.41782, 0, 0.212978, 0, 0.8178, 0,
0],
S1 [ 0, -0.772176, 0, 0, 0.0419199, 0,
0],
S2 [ 0, 0, -1.45948, 0, 0, 0,
0.7046],
S5 [ 0, 0, 0.438294, -0.559132, 0.404227, 0,
0],
S6 [ 0.41782, 0, 0.330554, 0, -1.26395, 0.616207,
0.346726],
S8 [ 0, 0.772176, 0, 0, 0, -0.616207,
-0.140416],
S9 [ 0, 0, 0.320719, 0.559132, 0, -0.616207,
-2.03405]]
|
|
I can confirm I am getting differences in the two jacobians with Matthias'
model.
|
|
@luciansmith Please use the model as test model. |
|
An observation, If I move two of the setvalue calls to the model itself (ie
in the antimony script):
vol_point = 0.1 # [liter]
f_fluid = 0.2 # [dimensionless]
f_fat = 0.3 # [dimensionless]
f_tissue = 1 - f_fluid - f_fat
Vext = vol_point * f_fluid
Vli = vol_point * (f_fat + f_tissue)
Both jacobians are the same and equal to the original reduced jacobian.
I'll continue looking.
|
|
I have computed the first entry in the jacobian manually and I get the same
answer that getFullJacobian and getReducedJacobian get. (Assuming I moved
Vext and Vli to the main program).
|
|
OK! This is now fixed with #912 The problem was that the compartment volumes weren't being synced with the species values when trying to deal with conserved moieties. If you didn't have to worry about that (i.e. the reduced form) there was no error. But now, both should match! Wheels available at https://dev.azure.com/TheRoadrunnerProject/roadrunner/_build/results?buildId=1103&view=artifacts&pathAsName=false&type=publishedArtifacts As a bonus, the manylinux build is now compiled in release mode, which should make it faster (though I don't think the llvm bit is affected). |
|
Thanks. Everything works. Thanks for the fast bugfix. |
|
Good to hear it.
|
Hi all,
I want to access the changes in reaction rates depending on species.
As I remember this is the information available via the Jacobian, i.e.
dv_i/dS_kis theJ_ik.So I tried to access the reduced Jabobian, but get a matrix which has species ids for rows and columns (and not rate ids for rows and species ids for columns)?
Is this correct? I would expect reaction ids in the rows here? What do the entries mean in this square matrix? And how would I access the changes of reaction rates with species?
When trying to access the full Jacobian, I get
so that I cannot access the values.
Code and model attached:
apap_core.zip
The text was updated successfully, but these errors were encountered: