Skip to content

Lab 1.1: Model formulation

The purpose of this lab is to learn how to formulate different types of models (ODEs and DAEs) based on previous information about the biological system, and how to analyze the different behaviours in the resulting model.

Throughout the lab you will encounter different "admonitions", which as a type of text box. A short legend of the different kinds of admonitions used can be found in the box below:

Different admonitions used (click me)

Background and Introduction

Useful information

Guiding for implementation

Here you are tasked to do something

Reflective questions

Before you start with this lab, you are expected to have:

  • Matlab running on your computer (or a computer in the computer hall)
  • A C-compiler installed in Matlab
  • Downloaded and installed IQM tools in Matlab
  • Downloaded and extracted the scripts for Lab 1.1: Lab1.1 files

If you haven't done this, go to Resources and follow the instructions before you start with this lab.

Create a mathematical model of fat metabolism

In systems biology, one common problem is to go from a scientific question about the mechanisms of a biological system, such as the fat metabolism, to formulating a mathematical model. Here, you are set to formulate, simulate, and analyze a mathematical model of the fat metabolism based on a hypothesis of how the metabolism works.

All files needed in this example are located in the folder fatMetabolism inside the Lab1A files.

First, we formulate the mathematical model.

Task 1: Formulate the mathematical model of fat metabolism
  • Read the model description
  • Formulate the model equations in a model file in matlab
Click to expand the hypothesis description

Figure1

Figure 1: Hypothesis of fat uptake from a meal. Created with BioRender.com

Figure 1 describes the hypothesis of fat uptake during a meal. We are interested in the dynamics of triglycerides (TAG) and free fatty acids (FFA). TAG consists of three FFA and one glycerol, and it is in this form that fat usually is stored in food and in the body. FFA is released when the body needs to utilize the energy from the fat - either to store it in another place or to use it as energy source.

All concentrations of meal-derived TAG and FFA in the body are 0 before the meal. The meal starts at time 0 and consists of 6000 \(\mu\)mol TAG.

The TAG from the meal is taken up in the blood with a rate of 0.015 min\(^{-1}\).

The TAG in the blood is then hydrolysed by the enzyme lipoprotein lipase (LPL), that cleaves the TAG into three FFA. The action of LPL is saturated by the concentration of TAG in the blood (\(TAGp\)) and can be described by Michaelis Menten kinetics, where the limiting rate \(Vmax = 17 \mu mol/L/min\) and the Michaelis menten constant \(km = 100 \mu mol/L\): \(vLPL = \frac{Vmax*TAGp}{km+TAGp} \mu mol/L/min\).

The FFA in plasma can be transported into the adipose tissue through diffusion. The diffusion rate pp is assumed to be 0.08 min\(^{-1}\). The amount of FFA that is transported from blood to adipose tissue and the other way around is dependent on the concentrations in tissue and blood, as well as the volumes in the tissue and blood respectively. The plasma volume is approximated to 3 L and the volume of adipose tissue is 1.2 L.

Finally, the FFA in the adipose tissue can be stored inside the adipose tissue, and this happens with a reaction rate around 1 min\(^{-1}\).

The concentration of both TAG and FFA in plasma after a meal are expected to increase quickly after the meal, to then decrease as the TAG is cleaved into FFA and the FFA is taken up and utilized in the tissue. The meal response is expected to last up to 6 hours (360 minutes) until the fat from the food is fully taken up in the tissues.

Implement model equations based on the interaction graph

When you have read the hypothesis and looked at the interaction graph, the next step is to translate the interaction graph into equations, in this case ordinary differential equations (ODE). Create a text file and write down the equations in the following IQM toolbox format. Save the file as fatModel.txt.

Equation format

The text file with the model equations (model file) should have the following format (see ExampleModel.txt): Note that the text file must contain all headers (beginning with “**** “)

********** MODEL NAME
ExampleModel
********** MODEL NOTES
********** MODEL STATES
d/dt(A)= r1 - r2                                                           
A(0) = 1
********** MODEL PARAMETERS
k1 = 1
k2 = 5
********** MODEL VARIABLES
variable1 = A*10
********** MODEL REACTIONS
r1 = k1
r2 = A*k2
********** MODEL FUNCTIONS
********** MODEL EVENTS
********** MODEL MATLAB FUNCTIONS

How to model food uptake

You can implement the food uptake as follows:

d/dt(TAGfood) = -foodUptake
TAGfood(0) = 0
foodUptake = TAGfood*kf

Where TAGfood is the amount of triglycerides that the food contains in \(\mu\)mol, \(foodUptake\) is the reaction of TAG transferred from the food into the blood, and \(kf\) is a parameter describing the rate of food uptake in min\(^{-1}\).

Note: Since the TAG in plasma is expressed as a concentration (\(\mu mol/L\)) and not just an amount in \(\mu mol\), the uptake of food into plasma is dependent on the volume of the plasma: \(foodUptake/Vplasma\).

How to model diffusion

One can think of diffusion as two reactions happening at the same time: one reaction \(diffPT\) transporting FFA from the plasma to the tissue, and one reaction \(diffTP\) transporting FFA from the tissue to the plasma. The diffusion has the same diffusion rate \(pp\) min\(^{-1}\) in both directions, but the amount of FFA being transported in each direction per minute depends on the amount of FFA in the respective place. The amount can be calculated from the concentration in \(\mu\)mol/L times the volume in \(L\).

The net diffusion (\(\mu\)mol/minute) in the direction from plasma to tissue can then be calculated as the difference between the transport to the tissue and the transport from the tissue: \(diffPT - diffTP\).

All these above calculations can be treated as variables in the model.

Finally, the diffusion reaction \(vdiffusionToTissue\) in \(\mu\)mol/L/minute can be calculated by taking into account the volume of the tissue: vdiffusionToTissue = (1/Vtissue)*netDiffusionToTissue. The diffusion reaction in the other direction, \(vdiffusionToPlasma\), can be calculated in a similar manner with respect to the volume in plasma: vdiffusionToPlasma = (1/Vplasma)*(-netDiffusionToTissue)

Note: Remember that in the conversion from TAG to FFA, the LPL reaction results in 3 FFA for every TAG.

Was the model you created correctly implemented? A first step to check this is to simulate the model and visually inspect if the simulation behaves as expected.

Task 2: Simulate the created model

To simulate and plot the model, we will use three scripts: a main script mainFat.m that will call your other scripts, a simulateExperiment.m that simulates the model, and plotSimulation.m that creates figures with plots of the model simulation. In this task, we will start by creating a simulation:

  • Compile the model in mainFat.m
  • Implement the simulation function simulateExperiment.m
  • Call the simulation function simulateExperiment.m from the main script mainFat.m
Compile the model

Open the script mainFat.m in Matlab. From this script, we will call all other functions needed to simulate and analyze the model.

To be able to simulate the model, we first need to convert your equations into a model object, so that you can then compile the model into a MEX-file. The MEX-file is essentially a way to run C-code inside MATLAB, which is much faster than ordinary MATLAB code.

In the mainFat.m script, add code to compile your model:

modelName='fatModel';
model=IQMmodel([modelName '.txt']);
IQMmakeMEXmodel(model) %this creates a MEX file which is ~100 times faster
[paramNames,params] = IQMparameters(modelName); %retrieve a vector with the parameter names and values from the model file (this is needed later)
Note: The compiled file will have a filename based on what is written under ****** MODEL NAME, and an ending based on what OS it is compiled on. In this case it should be fatModel.mex[w64/maci64/a64].

When you have compiled the model, it is time to do a simulation.

Implement the simulation function

To simulate the model, we use the function simulateExperiment.m. This simulation function should take (params, modelName, time, and IC) as inputs and returns the struct simulation as output in the following structure:

sequenceDiagram
MainFat ->> simulateExperiment: params, modelName, time, IC
simulateExperiment ->> MainFat: simulation

The blue boxes represent a function/script with the name written in it. The black arrows represent the direction a set of arguments is sent between two functions. In the above example, the script MainFat passes the arguments params, modelName, time, IC to the function simulateExperiment, that returns simulation back to MainFat.

In the simulation function ´simulateExperiment.m´, implement the actual simulation of the model: Since we have compiled the model, we can now simulate the model by directly calling the MEX file. In practice, we convert the name of the model in the form of a string (in modelName) into a callable function using str2func:

model = str2func(modelName); % Defines a callable function from the string value in modelName
sim = model(time, IC, params); % Simulate the model

Here, sim is a struct containing the resulting simulation, fatModel is the name of the model given when you compiled the model, time are the desired time points to simulate, IC is the initial conditions, and params are the values of the parameters in the model.

Note: it is also possible to skip IC and/or params as inputs. The model will then instead use the values written in the model file (when the MEX-file was compiled). In practice, this is done by passing along [] instead of e.g. IC.

Call the simulation function from the main script

Now, it is time to add the call to the simulation function in the mainFat script.

To simulate the food intake starting at time 0, we first set the initial condition of the state \(TAGfood\) to the amount of fat that the food contains (6000 \(\mu\)mol, see hypothesis description), while the other initial conditions are 0 since no food-based fat is assumed to be in the plasma or adipose tissue when the meal starts:

% Assumes that the state TAGfood is the first state in the modelfile.
IC = [6000 0 0 0]

Set time to a vector in minutes. The time should start at 0, and end at a time that is long enough time to see the meal response.

time = [0:endtime]; % you need to define a reasonable endtime in minutes

Then, we use the simulateExperiment function to get the simulation.

simulation = simulateExperiment(params, modelName, time, IC)

The resulting simulation should be a struct containing ie simulation.statevalues, which are the values of all states in the model for all simulated time points.

Using the simulation function, we can now plot and analyze the model simulations.

Task 3: Plot and analyze the model simulations

In this task, you should create the following plots:

  • Plot an overview of everything in the model
  • Plot the diffusion reactions
    • Set the reaction rate of the utilization of FFA in the adipose tissue to 0, and plot the diffusion reactions again
  • Plot the LPL reactions
    • Change the value of Vmax, and plot the LPL reactions again
    • Change the value of Km, and plot the LPL reactions again

Note: For each plot task, there is a prepared Matlab function available for you in the Lab1A files.

First, we will plot and analyze the general behaviour of the model.

Plot an overview of everything in the model

Call the script ´plotSimulation.m´ from your main script to plot all states, reactions, and variables in the model. ´plotSimulation.m´ takes (params, modelName, time, IC) as inputs but does not return any outputs:

sequenceDiagram
Main ->> plotSimulation: params, modelName, time, IC
plotSimulation ->> Main: [ ]

Look at the plot with all the model states. How does the meal response look? Is it reasonable compared to the hypothesis?

If the model is implemented correctly, the following questions should be answered with a yes:
  • Is the TAG in the food disappearing as it is taken up into the blood?
  • Does the TAG and FFA in the blood increase directly after the meal? The peak of TAG in plasma should be around 56 minutes at 440 \(\mu\)mol/L TAG, the peak of FFA in plasma around 75 minutes at 552 \(\mu\)mol/L, and the peak of FFA in tissue around 75 minutes at 102 \(\mu\)mol/L.
  • Is the fat utilized by the body? How long does it take until all fat is utilized?

Note: If you detect any errors in your model formulation and make changes in the model file, you need to re-compile the mex file. It might be a good habit to always compile the mex file before you do any simulations.

Look at the plot with all the model reactions.
  • Identify the plots of the reactions that correspond to the mass action kinetics. Which are they?
  • Identify the reactions that correspond to the diffusion between blood and adipose tissue. Which are they?
  • Identify the reaction that correspond to the Michealis-Menten kinetics. Which is it?

To further analyze the diffusion and LPL reactions, we will specifically plot their reaction rates.

Plot and analyze the diffusion reactions

Call the function plotDiffusion.m from your main script to plot the reactions connected to the diffusion. ´plotDiffusion.m´ takes (params, modelName, time, IC, netDiffusionToTissue,vdiffusionToTissue, vdiffusionToPlasma, diffPT, diffTP) as inputs but does not return any outputs:

sequenceDiagram
Main ->> plotDiffusion: params, modelName, time, IC, netDiffusionToTissue,vdiffusionToTissue, vdiffusionToPlasma,diffPT,diffTP
plotDiffusion ->> Main: [ ]

The names netDiffusionToTissue,vdiffusionToTissue, vdiffusionToPlasma, diffPT, and diffTP should be strings (for example, the function input for netDiffusionToTissue can be defined as 'netDiffusionToTissue').

Make sure that you provide the same names of the reactions as in your model file as inputs to plotdiffusion:

  • netDiffusionToTissue: net amount of FFA diffusion in the direction towards the tissue (ie diffPT-diffTP).
  • vdiffusionToTissue: reaction rate of the net diffusion to the tissue (\(\mu\)mol/min/L)
  • vdiffusionToPlasma: reaction rate of the net diffusion to the plasma (\(\mu\)mol/min/L)
  • diffPT: amount of fat from plasma to tissue (\(\mu\)mol/min)
  • diffTP: amount of fat from tissue to plasma (\(\mu\)mol/min)
Look at the diffusion plots.
  • Which way does the net diffusion go to? Does this make sense?
  • Why is the reaction rates proportionally different from the amount of TAG that is being transported?
  • Does the diffusion ever stop? Why?
To test if the diffusion is correctly implemented, set the reaction rate of the utilization of FFA in the adipose tissue to 0.
% assumes that the parameter governing the fat utilization is named ku
params(strcmp(paramNames,'ku')) = 0; %turn off fat utilization

Now, the FFA will not disappear from the system. Instead, all FFA from the meal will accumulate in the adipose tissue and in the blood. Eventually, the model should reach a steady state where the net diffusion between tissue and plasma is 0. Run the mainFat script again, and look att the diffusion plots.

Look at the diffusion plots after setting the utilization to zero
  • Does the diffusion reach a steady state after the meal, when the net diffusion is 0? If not, go back to your model formulation of the diffusion reactions! (why?)
  • In the steady state, is the amount of FFA diffusing into the tissue the same as the amount of FFA diffusion out of the tissue? If not, go back to your model formulation of the diffusion reactions! (why?)
  • In the steady state, is the reaction rate of FFA diffusing into the tissue the same as the reaction rate of FFA diffusion from tissue to plasma? If yes, go back to your model formulation of the diffusion reactions! (why?)

Note: Before continuing, remember to set the reaction rate of the utilization of FFA in the adipose tissue back to 1 again!

Plot and analyze the LPL reactions

Call the function ´plotLPLkinetics.m´ from your mainFat script to plot the reactions connected to the diffusion. plotLPLkinetics.m takes (params, modelName, time, IC, plasmaTAGname, LPLreactionname) as inputs but does not return any outputs:

sequenceDiagram
Main ->> plotLPLkinetics: params, modelName, time, IC, plasmaTAGname, LPLreactionname
plotLPLkinetics ->> Main: [ ]

The state- and reaction names (plasmaTAGname and LPLreactionname) should be provided as strings.

Make sure that you provide the same names of the state TAG in plasma (plasmaTAGname) and reaction of LPL (LPLreactionname) as in your model file as inputs to plotLPLkinetics.

Look at the plot of the velocity of the LPL reaction vs the concentration of the substrate TAG in plasma.
  • Can you find the substrate concentration where the velocity is half of Vmax? What parameter does that velocity correspond to?
  • Does the rate of the LPL reaction ever reach Vmax?
  • Compare to the plot of the LPL reaction vs time and the plot of the concentration of TAG in plasma vs time - can you find the corresponding time point where the velcity is Vmax/2?
  • Change the value of the km parameter (ie to 10 or 1000 instead of 100) and re-run the main script. How does that affect the plot of LPL kinetics?
    % assumes that the km parameter is named km in the model file
    params(strcmp(paramNames,'km')) = 10;
    
  • Change the value of Vmax (ie to 30 or 10 instead of 17) and re-run the main script. How does that affect the plot of LPL kinetics?

To see how the model formulation affects the model simulation, we can change the kinetics of the LPL reaction and se how that affects the reaction rate.

Task 4: Change the kinetics of the LPL reaction

Sometimes one needs to try different types of kinetics for the reactions to see which kinetics that is needed to be able to describe the biological system and agree with data. Now, change the kinetics for the LPL reaction to see how the dynamics of the meal response changes.

Change the dynamics to the following kinetics and plot the results:

  • Mass action kinetics
  • Hill kinetics, with different values for n
Different types of kinetics

There are several types of kinetics that can occur in biological systems. Here, we will cover thee of the most common examples:

Mass action kinetics is the most basic type of kinetics, where the reaction rate \(v\) is directly proportional to the concentration of the reactant \(A\): \(v = k*[A]\)

Michelis Menten kinetics describes a reaction rate \(v\) that is saturated when there is a lot of substrate, and was originally used to describe enzyme kinetics. The Michaelis Menten kinetics can be expressed as: \(v = \frac{Vmax*[A]}{k_{m}+[A]}\) where the reaction rate \(v\) is saturated at the limiting rate \(Vmax\) if there is enough substrate \(A\). The Michaelis Menten constant \(k_{m}\) corresponds to the substrate concentration where the reaction rate is half of \(Vmax\).

Hill kinetics describes a sigmoidal response and was originally used to describe the binding of oxygen to hemoglobin in the blood. Hill kinetics can be used to describe receptor-ligand or enzyme-substrate kinetics in biochemistry or pharmacology. The Hill equation can be expressed as: \(v = \frac{Vmax*[A]^{n}}{k_{a}^{n}+[A]^{n}}\) where \(k_{a}\) is the substrate concentration where \(v = Vmax/2\) and the inflection point of the sigmoid response is found. The Hill coefficient \(n\) governs the sensitivity to the amount of substrate A: if \(n\) > 1 the binding of one substrate facilitates binding of more substrates, but if \(n\) < 1 the binding of one substrate reduces the affinity for other substrates. What happens when \(n\)=1?

Change the kinetics of LPL to mass action kinetics

Create a copy of your model file fatModel.txt, and name it fatModelMassAction.txt. In your new model file, update the model name and change the reaction vLPL from Michaelis Menten kinetics to mass action kinetics.

In your mainFat.m script, change the modelName to your new model fatModelMassAction, and run the script.

Now, we are only interested in the plot of the LPL kinetics. If you want, you can comment out the two lines calling the plot functions plotSimulation and plotDiffusion by adding a % before each line.

How does the meal dynamics change?

Look at the plot of the LPL kinetics. What is the difference between Michaelis Menten and mass action kinetics?

Change the kinetics of LPL to Hill kinetics

Create a copy of your original model file fatModel.txt, and name it fatModelHill. In your new model file, update the model name and change the reaction vLPL from Michaelis Menten kinetics to Hill kinetics by adding the Hill coefficient \(n\): \(vLPL = \frac{Vmax*TAGp^{n}}{km^{n}+TAGp^{n}}\)

First, set n = 3. Then, try at least two other values for \(n\): one where \(n<1\) and one where \(n=1\).

In your mainFat.m script, change the modelName to your new model fatModelHill, and run the script.

How does the meal dynamics change?

Look at the plot of the LPL kinetics.

  • What is the difference between Hill, Michaelis Menten, and mass action kinetics?
  • How does the value of the Hill coefficient \(n\) affect the kinetics?

Create a mathematical model of the cardiovascular system

Not all systems can be described using ordinary differential equations - there are several other types of equations that might be needed for different systems. In this example, to model the cardiovascular system, we need to add algebraic equations with time-varying variables without their derivatives. Thus, in this example we will use differential algebraic equations (DAE) instead of ODEs.

Differential algebraic equations (DAE)

DAEs are equations where ordinary differential equations are not enough to describe the system. Instead, additional algebraic equations are needed, which contain other time-dependent variables without their derivative: \(F(x'(t),x(t),t) = 0, x'(0) = x'_{0}, x(0) = x_{0}\), where \(x\) not only contains the state of the model but also other variables.

DAEs appear in many situations, such as electrical circuits, chemical processes, or mechanical systems like the motion of a pendulum.

The initial conditions of both \(x(t)\) and \(x'(t)\) can be set in DAEs. However, the initial conditions always need to be consistent with the equations you are solving, which can be especially tricky in DAEs. There are helper functions to find consistent initial conditions as close as your wanted ones as possible, and some solvers do this automatically.

Note: Some DAEs (ie with low differential index*) can easily be re-written as ODEs and solved with ODE methods - but not all.

*Differential index: The minimum number of times the DAE has to be differentiated to get an ODE. A semi-explicit DAE with index 1 could be defined as \(x'(t) = f(x(t),y(t),t), 0 = g(x(t),y(t),t)\). Be aware of that if a DAE is differentiated, the constrains might disappear.

All files needed in this example are located in the folder cardiovascularSystem inside the Lab1A files.

Cardiovascular models

Electrical circuits are often used as an analogy when creating models of the cardiovascular system, where the cardiovascular fluid dynamics are represented by electrical components: the blood flow is represented by a current I, pressure differences by a voltage V, vessel resistance by a electrical resistance, compliance of elastic vessels by conductance, and the inertia of the blood by inductance. Thus, by drawing an interaction graph of a simple circuit representing a part of the cardiovascular system, one can use electrical and physical laws such as Kirchhoff's current and voltage laws and Ohm's law to derive mathematical equations of the cardiovascular system.

Cardiovascular fluid dynamics Electrical circuit
Figure2 Figure2
Vessel resistance Electrical resistance
Blood flow Current
Pressure difference Voltage difference
Compliance (distensibility) of elastic vessels Capacitance
Blood inertia Inductance

Figure2

Figure 2. An example of a cardiovascular model from Casas et al 2018, where each part of the heart to the left corresponds to a part in the electrical circuit to the right.

Based on a description of the cardiovascular model, we will formulate the model equations.

Task 5: Implement the model equations of a simple model of the cardiovascular system

In the file cardiovascularModel.m, write down the model equations, model parameters, and model input in the indicated places.

Since the IQMtools that you used in the previous example does not include a DAE solver, we will instead formulate the model directly in Matlab.

Click to expand the model description

Figure2

Figure 2. The cardiovascular model, where the cardiovascular system with the heart and elastic aorta (top) is represented by an electrical circuit (bottom).

This simple model represents the blood flow and blood pressure in the aorta, where in general:

Q: Blood flow (mL/s)

R: Resistance (mmHg·s/mL), representing the resistance in blood vessels

L: Inertance (mmHg*s\(^{2}\)/mL), representing the inertia of blood flow

C: Compliance (mL/mmHg), representing the windkessel effect of the aorta which expands and stores blood during systole, and relaxes during diastole

Retrieve the equations from the circuit

To retrieve the equations from the circuit, we use electrical and physical laws:

1. Definition of Ohm's law, inductance, and compliance

Ohm's law: \(\varDelta U = R*I <-> \varDelta P = R*Q\)

The pressure over the inductance is given as: \(\varDelta P = L * \frac{dQ}{dt}\).

Compliance defines the stiffness of the cardiovascular system - how much volume change is generated by a change in pressure: \(C = \frac {\varDelta V}{\varDelta P}\). The rate of change in volume over a segment is the blood flow Q ml/s, which gives us:

Eq 1. \(Qdiff = Ca*\frac{dPa}{dt}\)

where \(Pa\) is the pressure over the compliance \(Ca\).

2. Kirchhoff's voltage law

All voltages (or blood pressure differences) in a loop should sum up to 0, which in combination with Ohm's Law and the definition of inductance gives:

Eq 2. \(Paorta - Rao*Qao - Lao* \frac{dQao}{dt} = 0\)

Where the pressure in the aorta is:

Eq 3. \(Paorta = Pa + Ra*Qdiff + La* \frac{dQdiff}{dt}\)

3. Kirchhoff's current law

All currents (and blood flow) in a node should sum up to 0:

Eq 4. \(Qheart - Qdiff - Qao = 0\)

Look at the equations
  • Which are the time-varying variables in the model that we will solve for?
  • Is it possible to re-write the DAEs to ODEs, with one equation defining each derivative, in a simple way?
Read the model hypothesis above, and write down the model equations in cardiovascularModel.m

Open the model file cardiovascularModel.m and write down the four model equations from the model description, the model input, and the model parameters in the indicated places.

Equation format

To solve DAEs, we cannot use the IQM toolbox. Instead, we will use the matlab solver ode15i and define the equations directly in Matlab.

For the example equation system (the Robertson problem), the model formulation is:

\(y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)\)

\(y(2)' = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2\)

\(y(3)' = 3e7*y(2)^2\)

We can reformulate the equations as DAE:

\(0 = y(1)' + 0.04*y(1) - 1e4*y(2)*y(3)\)

\(0 = y(2)' - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2\)

\(0 = y(1) + y(2) + y(3) - 1\)

And finally define the three equations in matlab as:

equations = [ yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)
yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2
y(1) + y(2) + y(3) - 1 ];

Where each row corresponds to one equation, y are the time-dependent variables, and yp are the first derivatives of the time-dependent variables.

Note that the Robertson problem is just an example - you task is to implement the four equations of the cardiovascular model described in the model description above.

For our cardiovascular model, we can define y as

Paa = y(1);
Qao = y(2);
Paorta = y(3);
Qdiff = y(4);
and we can define yp as
dPaa = yp(1);
dQao = yp(2);
dPaorta = yp(3);
dQdiff = yp(4);

For the plot script to work, \(Paorta\) needs to be y(3)

Parameters in the model

The following values for the parameters in the model represent normal values in the cardiovascular system:

Ca =  1.4; 
La = 0.0001;
Ra = 0.1;
Rao = 1.2;
Lao = 0.0005;

Input Qheart to the model

The input to the model is the blood flow from heart Qheart, which we model as a sinus curve:

Qheart = ceil(max(0,k_syst-mod(t,T)/T))  * 300.*sin((pi*(mod(t,T)/T))/k_syst);

Where k_syst = 0.4.

In diastole, the heart relaxes and the aortic valve is closed, which means that the blood flow from the heart is 0. This is modelled by the first part of the expression, which is 0 when the valve is closed and 1 when the valve is open. In systole, the aortic valve is open and blood is pumped out, which we model as a sinus curve with the amplitude 300 mL/s that depends on the ratio k_syst of the length of systole (tsyst) to the length of the full cardiac cycle (T): k_syst = tsyst/T.

The cardiac cycle is simulated to start at the start of systole at t=0.

What is the difference between an ODE and a DAE?
  • What makes this model a DAE instead of ODEs?

When you have written down all the equations in the model file, it is time to simulate the model and plot the results.

Task 6: Simulate, plot, and analyze the cardiovascular model

Run the existing mainCardiovascular.mscript to simulate tree heartbeats with the model and to plot the resulting simulated blood pressure.

Look at the resulting plot
  • Blood pressure is usually measured as the systolic (the highest value, during systole) and diastolic (the lowest value, during diastole) pressure. What is the systolic and diastolic blood pressure simulated by the model? If the equations are correctly implemented, the systolic pressure should be 118 mmHg and the diastolic pressure should be 71 mmHg.

When aging, the blood vessels get stiffer. The stiffening results in a lower compliance of the blood vessels. Do a model prediction of this process by reducing the value of the model parameter Ca and running the mainCardiovascular again.

Look at the resulting plot
  • In what way did the systolic and diastolic pressure change with the siffening of the vessels?

Questions for lab 1.1 - model formulation:

In your lab report, you should answer the following questions:
  • Show one of your plots of Michaelis-Menten kinetics. Describe how the parameters km and Vmax govern the appearance of the plot.
  • Show your plots of the diffusion with and without fat utilization. What is the difference? Explain how you can tell that the diffusion equations are correctly implemented in your model.
  • What is the difference between an ODE and a DAE? Use the equations from the two example models from the lab and describe why they are ODEs or DAEs respectively.