Modelling a system of ODEs and need to add inputs after t=0?

Hi all,

I'm testing a pharmacokinetic model in MATLAB, it's very stiff so I'm currently using ode15s, the system is comprised of 7 interconnected ODEs representing various 'compartments' within the body.

My issue is as follows: I'm attempting to forward model the system based on a fitted parameter set, in order to model is accurately I need to mimic the inputs as done in the original in vivo trials. There in one input at t=10s and two further inputs at t=60s.

I've attempted a work around, just attempting to brute force it in the function but that led the solver to complain about infinitely small step sizes so I assume I broke something.

Any suggestions would be fantastic! Thank you.

Answer 1

[Star Strider]

I am not exactly certain what you are doing or what your pharmacokinetic model is. Not knowing more, I would do something like this (assuming ‘t’ is the time vector and the integrated values are stored in matrix ‘y’):

Use first set of initial conditions ‘ic’

tspan = linspace(0, 10, 10);

Run ode15s with your ODE function

ic = y(end,:);
t0 = t(end);
ic(appropriate_compartment) = dose(1);

Use new set of initial conditions ‘ic’

tspan = linspace(t0, t0+60, 60);

Run ode15s with your ODE function

ic = y(end,:);
t0 = t(end);
ic(appropriate_compartment) = dose(2);

... continuing for as many steps as necessary. Adjust the number of elements in ‘tspan’ with the third argument of linspace. You will need to vertically concatenate the ‘t’ vectors and ‘y’ matrices from each step (loop iteration) to create the complete simulation.

I’m simply sketching this, corresponding to similar problems I’ve worked on.

Apparently the SimBiology has a set of functions that will make alll this straightforward, although I have no experience with it. (I don’t do such simulations frequently enought to justify buying it.)

[Arthur Goldsipe]

As Star Strider mentioned, SimBiology is designed to make it easy and fast to simulate exactly this sort of model. If that's an option you want to learn more about, I suggest creating another MATLAB Answers question asking specifically about building PK models in SimBiology.

[Daniel Bending]

Hi all, thank you for the replies! SimBiology does sound interesting, though at the moment I'm just trying to verify a model made in monolix (It worked in August, doesn't work now!).

I think what you're suggesting is effectively a superposition of two seperates instances of the model, right? I did think about that, but I believe the concept behind solving stiff models is a time step that is variable? I suppose I can force a very small time step so they variable develope correctly..!

I'm not entirely sure that option will work either, as there are a number of varaibles that are calibrated around basal levels (i.e. the level of glucose etc) in the blood pre-experiment. As a result if I attempted to run the simulation with different initial conditions I would need the basal level input into each instance which, when added together would have twice the total amount. In theory I could just subtract the extra basal level but as there are non-linear terms and control terms (i.e. the amount of x in the blood causes y to happen without material transfer) this might not work anyway!

Thanks again!

[Star Strider]

My pleasure.

I’m not certain what your model is or the reservations you might have about the procedure I outlined. I’m familiar with the physiological adaptations that naturally occur with varying plasma glucose concentrations (including counter-regulatory hormones and receptor changes), and a sufficiently nonlinear model would of course take those into consideration and model them appropriately. I have not done such modeling because estimating changing parameters is extremely difficult. From a control theory perspective, a model must be both observable and controllable.

I have no idea what your model is or what it entails, onlly that it apparently has a set of identifiable (assumed non-changing) parameters. The procedure I’ve outlined is the typical way a relatively static set of ODEs would be simulated with different inputs.

I’ve posted a number of solutions here on kinetic model parameter estimation and simulation using various MATLAB ODE solvers.

[Daniel Bending]

Thank you again for your very prompt reply.

In regards to the model, I've recently been given the model, the population data it's been fit to and the best parameter fit from the time. Unfortunately it completely fails to run in the software it was made in so I'm attempting to test it in MATLAB.

Actually stratch everything I previously said, what you've suggested is genius and I wasn't reading it properly! You run the simulation until t=(time of first input) take the output of this section and put it back as the new initial conditions plus the inputs! That's fanstic! I was completely misunderstanding. Yes that would absolutely work! Thank you!

[Star Strider]

As always, my pleasure!

If you care to post the differential equations and parameters, I will do my best to get it to work with ode15s. With dosages and the input compartment, I might be able to make it work with those as well. (I can’t promise anything because I’ve not seen the differential equations or the data that was used to estimate them.)

[Daniel Bending]

Aha, you are a star! Alas, I've probably said too much already, it's for an industry sponsored masters research project. I've also managed to get it to work in the time, your advice was very clear once I'd opened my mind enough to actually read it!

It's clear now that the paramatization is completely off, there's an input of 11 nmol/L of insulin and it instantly shoots up into the millions, I think not..!

[Star Strider]

Thank you!

I understand. If at some point you’d like my help in solving it, post back here with the model and the data and I’ll do my best to identify it correctly and simulate it.

In the interim, see: Parameter Estimation for a System of Differential Equations. That’s likely the most concise of such posts. There are others, including some that use the ga (genetic algorithm) to estimate the parameters.

[Sietse Braakman]

Hi Daniel,

One of the problems you are running into is the fact that you need to restart the ODE solver, every time you dose/perturb the model. SimBiology can really help here, as it takes care of the ODE solver restarts and makes dosing much easier. The function sbiofit (or a fit program in the Model Analyzer App) will allow you to set up parameter estimations easily and iterate through different algorithms, try out different error models. Primarily, it constructs the objective function for you.

PM me if this is of interest - I can help you get started with SimBiology.

Best,

Sietse