Resultados de
Automating Parameter Identifiability Analysis in SimBiology
Is it possible to develop a MATLAB Live Script that automates a series of SimBiology model fits to obtain likelihood profiles? The goal is to fit a kinetic model to experimental data while systematically fixing the value of one kinetic constant (e.g., k1) and leaving the others unrestricted.
The script would perform the following:
Use a pre-configured SimBiology project where the best fit to the experimental data has already been established (including dependent/independent variables, covariates, the error model, and optimization settings).
Iterate over a defined sequence of fixed values for a chosen parameter.
For each fixed value, run the estimation to optimize the remaining parameters.
Record the resulting Sum of Squared Errors (SSE) for each run.
The final output would be a likelihood profile—a plot of SSE versus the fixed parameter value (e.g., k1)—to assess the practical identifiability of each model parameter.
I want to observe the time (Tmax) to reach maximum drug concentration (Cmax) in my model. I have set up the OBSERVABLES as follows (figure1): Cmax = max(Blood.lL15); Tmax_LT = time(Conc_lL15_LT_nm == max(Conc_lL15_LT_nm)); Tmax_Tm = time(Conc_lL15_Tumor_nm == max(Conc_lL15_Tumor_nm)); After running the Sobol indices program for global sensitivity analysis, with inputs being some parameters and their ranges, the output for Cmax works, but there are some prompts, as shown in figure2. Additionally, when outputting Tmax, the program does not run successfully and reports some errors, as shown in figure2. How can I resolve the errors when outputting Tmax?
Hi All,
I'm currently verifying a global sensitivity analysis done in SimBiology and I'm a touch confused. This analysis was run with every parameter and compartment volume in the model. To my understanding the fraction of unexplained variance is 1 - the sum of the first order variances, therefore if the model dynamics are dominated by interparameter effects you might see a higher fraction of unexplained variance. In this analysis however, as the attached figure shows (with input at t=20 minutes), the most sensitive four parameters seem to sum, in first order sensitivities to roughly one at each time point and the total order sensitivies appear nearly identical. So how is the fraction of unexplained variance near one?
Thank you for your help!
Hi to everyone!
To simplify the explanation and the problem, I simulated the kinetics of an irreversible first-order reaction, A -> B. I implemented it in two independent compartments, R and P. I simulated the effect of a dilution in R by doubling at t= 0,1 the R volume. I programmed in P that, at t = 0.1, the instantaneous concentration of A and B would be reduced by half. I am sending an attach with the implementation of these simulations in the Simbiology interface.
When the simulations of the two compartments are plotted, it can be seen that the responses are not equal. That is, from t = 0.1 s, the reaction follow an exponential function in R with half of the initial amplitude and half of the initial value of k1. That is, the relaxation time is doubled. Meanwhile, in P, from t = 0.1, the reaction follows exponential kinetics with half the amplitude value but maintaining the initial value of k = 10. Without a doubt, the correct simulation is the latter (compartment P) where only the effect is observed in the amplitude and not in the relaxation time. Could you tell me what the error is that makes these kinetics that should be equal not be?
Thank you in advance!
Luis B.
Hi All,
I've been producing a QSP model of glucose homeostasis for a while now for my PhD project, recently I've been able to expand it to larger time series, i.e. 2 days of data rather than a singular injection or a singular meal. My problem is as follows: If I put 75g of glucose into my stomach glucose species any later than (exactly) 8.5 hours I get an integration tolerance error. Curiosly, I can put 25g of glucose in at any time up to 15.9 hours, then any later an error. I have disabled all connections to my glucose absorption chain, i.e. stomach -> duodenum -> jenenum -> ileum -> removal, to isolate the cause of this. I had initially thought it may be because I mechanistically model liver glycogen and that does deplete over time, but I've tested enough to show that that does nothing. My next test is to isolate the glucose absorption chain into a seperate model and see if the issue persists but I'm completely baffled!
These are the equations, to my eye there's no reason why there would be such a sharp glucose quantity/time dependence, they all begin at a value of 0:
d(Gs)/dt = -(kw*(1-Gd^14/(Igd^14+Gd^14))*Gs) #Stomach glucose
d(Gd)/dt = (kw*(1-Gd^14/(Igd^14+Gd^14))*Gs) - (kdj*Gd) #Duodenal Glucose
d(Gj)/dt = (kdj*Gd) - (kji*Gj) #Jejunal Glucose
d(Gi)/dt = (kji*Gj) - (kic*Gi) #Ileal Glucose
(The sigmoidicity of gastric emptying slowing term (^14) was parameterised off of paracetamol absorption data and appears to be correct!)
Thank you for your help, best regards,
Dan
Pre-Edit: I changed the run time to 30 hours and now I can't use the 75g input any later than 7.9 hours not 8.5 hours anymore!
Edit: This is how it appears at all times prior to it failing for 75g:
Join us for a hands-on workshop focussing on Global Sensitivity Analysis (GSA) for QSP using SimBiology. This complimentary session aims to give you theoretical background and hands-on experience with GSA methods.
🗓️ Event Details: June 25, 1-5pm CEST @PAGE 2024, Rome
🔖 Registration: Free, but space is limited. Early registration is recommended to ensure your place in this workshop.
For more details and to register, please visit https://www.mathworks.com/company/events/seminars/global-sensitivity-analysis-with-simbiology-4361208.html
Hi All,
I'm trying to produce some nice figures from the graphical user interface and have a set of local sensitivity analyses that I'd like to combine.
I have two inputs that vary the sensitivies of my system and would like to plot them on top of each other like:
Is there either a way to do this natively in simbiology (when you try and use 'keep results from each run' it plots them both as a time series) or to export the sensitivity data to the normal matlab programatic UI where I can combine them by hand?
Thank you for your help!
Hi All,
I'm trying to plot my observables in the analyzer and running into far more problems than expected. I have two species in my model that are invovled, blood.Insulin (pM) and blood.Glucose (mM), all I want to do is plot the ratio of these two (blood.Insulin/blood.Glucose (dimensionless)) along with my other species in Model Simulation, to compare it to the same ratio from my data.
First, there doesn't seem to be a way to directly add an observable to the logged states in 'Simulate Model', so I've tried to used 'Calculate Observable' based on the data from my last run (IVGTT.LastRun.results) but it says that units are required when unit conversion is enabled, but it should be dimensionless!
My next idea would be to make a non-constant parameter with a repeated assignment, but I feel like I should be able to do this without resorting to that?
Any help or ideas would be appreciated. Thank you, best regards,
Dan
I based my model construction on this PBPK model: PBPK by Armin Sepp. While this is a very convenient script for building a PBPK two-pore model, it's very incovenient for my application to have the species Units defined in molarity. Is there a convenient way to organically switch this model from molarity to grams (or any weight unit)?
Hi! I'm new to pk modeling and Matlab. Can someone guide me through how to conduct population pk modeling based on pk parameters from non-human primate studies? Much thanks!!!!
The title is resonably non-descript, but I can explain it easily:
Say I have an initial Emax model:
v = emax1*[G]^n1/(ec501^n1+[G]^n1)
And I want to place v inside of a second Emax model:
y = emax2*v^n2/(ex502^n2+v^n2)
Currently, I have the full function of v inside y, twice, it's very long and whilst I only need to get it correct once, for readability in the future I'd rather have it in form #2. I've played around with non-constant parameters but I need the steady state to be v, not the rate rule, and I haven't worked out how to make a parameter shift to a form like v, as an observation might.
Are there any recommended solutions or do I simply need to keep with having v fully expressed in y?
Thank you,
Dan
Hi All,
I'm currently attempting to implement a Hodgkin-Huxley-type model of membrane potential, ideally I would like a species that represents the membrane potential as its own distinct entity, so as the reference elsewhere. I've currently established a molarity-based work around but it would be great if I could set the units for the species as millivolt, but that throws an error.
Is there an established way to do this? I imagine I'm not the first person to be trying to model a voltage-gated ion channel!
Thank you for your help.
Hello,
I've looked around and I haven't found anything obvious about this, but is it possible to link to species/reactions, graphically, in a non-mass transfer sense? I have areas in my model where it would conceptually make sense to be able to see that species or reactions are linked, but if I link them in the standard way it demands that it be involved in the stoichiometry.
Perhaps some kind of dotted line, or similar?
Thank you, best regards,
Dan
Hello all,
I've been trying to shift my workflow more towards simbiology, it has a lot of very interesting features and it makes sense to try and do everything in one place if it works well..! Part of my hesitancy into this was some bad experiences handling units in the past, though this was almost certainly all out of my own ignorance, relatedly:
Getting onto my question.
In this model I have a species traveling around the body via blow flow, think a basic PBPK model. My species are picomolarities, if everything is already in concentrations, why is it necessary to initially divide by the compartment volume? i.e. 1/Pancreas below.
If my model dealt in molar quantities this would make a lot of sense, the division would represent the transition to concentrations. This, however, now necessitates my parameters be in units of liter/minute, which is actually correct, but I'd like clarification on why it's correct, ha!
Perhaps this is more of a modelling question than a simbiology question, but if there are answers I'd love to hear them. Thanks!
Hi All,
I'm attempting to put a set of simbiology global sensitivity analysis plots into my thesis and I'm running into some issues with the GSA plots. Firstly, the figures are very large, it would be quite beneficial to grab a set of the plots and arrange them myself, is there any documentation on how to mess around with the '1x1 Sobol' produced by sbiosobol? Or just GSA plots in general.
The second problem is that the results appear to be relative to the most sensitive parameter in that run. Is it recommended to have a resonably sensitive 'baseline' parameter in each run? I find it difficult to compare plots when a not so sensitive parameter is being recorded as near '1' for the whole run because it's being stacked against a set of very insensitive parameters. I.e. if i have multiple sets of GSAs due to a large model, how can I easily compare results? If I could do some single run through with every parameter that would be the ideal, I imagine, but then the default plot would be half a mile off the bottom of my screen, haha! Perhaps there is a solution to the first question that might help there?
Thank you for your help,
Dan
Hi all,
I've translated a model from another piece of software (monolix) into simbio programmatically to make use of your very easy global sensitivity analysis system.
It looks a little something like this, for a 'single' line example:
r1 = addreaction(model,'InsI -> InsP');
r1.ReactionRate = 'InsI*kip/vi'; %- is + panc
k1 = addkineticlaw(r1, 'Unknown');
Multiplied about 20 fold, as you can see I have included my volumes within the reaction rates myself (vi). The model functions perfectly and I have corrected the outputs at the end:
[time, x, names] = sbiosimulate(model,csObj,dObj1);
x(:,1) = x(:,1)/vi;
So that they are in concentration, as needed. However, when it comes to sensitivity analysis because I have corrected them post-model it is technically incorrect, it is analysing the absolute quantities. This is quite noticible in the sensitivity to the volumes.
Is there an easy fix to this, I've had to fight dimensionality with units in the past using simbio and I'd be great if there was some way of dividing a compartment output by a volume, for example. It is a functionality that exists in monolix, so I was hopeful it might here!
Thank you for your time.
EDIT:
I think I've worked it out, I had to refactor my model to operate in concentrations, just refitting it now. Now I should just be able to use unitless compartments.
I am trying to simulate model of blood lymphocyte count from a paper using Simbiology. The rate of in or out following circadian rythym is kp(t)=km + kb cos [(t-tpeak)*2pi/24] Where and how do I write the expression ? I dont think I can write in repeated assignment ?
I am currently facing a compatibility issue when attempting to load a SimBiology model created in MATLAB 2021a into MATLAB 2023a. Specifically, the code capture functionality does not seem to be working.
If anyone has encountered a similar situation or has insights on how to capture the code for a SimBiology model created in an older version of MATLAB, I would greatly appreciate your guidance. Are there any alternative methods or specific steps that can be followed to ensure successful code capture?
Thank you in advance for your time and expertise.
Dear community,
I would like to develop a fermentation model with 4 ODEs, one of which contains variable y. A "repeated assignment", e. g. y=5x+5, contains variable x that has been measured each second. These data (columns with time and corresponding value x in each row) are recorded in the Excel file.
Does anyone have any suggestion how to implement this in symbiology?
Thank you very much in advance,
Tetiana
Hi,
I want to develop a PK model based on some PK data. The PK data seems to display 2 peaks when one initial dose is given.
I would like to give one dose. A fraction of this dose (fr) is absorbed following the linear relationship - ka*Drug the other fraction (1-fr) is absorbed following a linear absorption (ka1*Drug) with a Tlag (it maybe a zero order). The fraction fr is unknown so it must be estimate.
Does anyone have can provide any suggestion to implement this in symbiology or provide a link where I can look?
Thank you very much in advance,
Ferran
