Writing the s-function as m-file

11:43 am

Hey!! Good morning! Back after indeed a very very long time! Had my 3rd semester exams; and then a long vacation of an entire week. Pampered myself for a while. Now back to the project. Pressure needs to build up for me to get back on track. Spoke to my cousin sister yesterday, and finally decided that I should start.

So, I have sort of implemented all that I have written in the previous post, only the most important one, the last block of pde solution isn't done. I tried to do it.

Firstly, followed the same steps they have given in the pdetool, parabolic part in the manual. It's missing the plotting code. So realised I need to understand the bloody code. If I don't understand it, it's senseless. To understand the code, .. I HAVE to blog! Else my mind wanders around. So, here we go..

First line in the code is:

[p,e,t]=initmesh('circleg'); % Create initial triangular mesh

As we all know, to apply the pde on every point of the geometry, we first divide the geometry into small triangles. This we call meshing. We apply the pde on every triangle and then integrate the result. This is the method of finite element I guess. Need to search a bit on what it exactly is, but I think that's it. And this function 'initmesh' initializes the creation of this mesh.
The function obviously needs the geometry as the input parameter, which we can define in different ways according to the format and flexibility of the function. The output of the function is the meshed data; p => point matrix with x and y coordinates of mesh points as the rows of the matrix, e => edge matric with starting and ending point indices, starting and ending point parameter values, edge number and left and right side subdomain number, t => triangle matrix with three corner points in anticlockwise direction and subdomain number.

Next,

[p,e,t]=refinemesh('circleg',p,e,t);  % Refine triangular mesh

Refinement usually bisects every edge of the specified triangles, thus creating 4 new smaller triangles out of every specified triangle. Hence the input to the function is the geometry under consideration, and the specified initial mesh and the output of the function is the p, e, t matrices of the refined mesh for that geometry.

Next,

u0=zeros(size(p,2),1);

size(p,2) returns the size of the 2-D p matrix and zeros(size(p,2),1) returns a 2-D zero matrix of the same size. We name it u0.

Next,

ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);  % Find indices and values of nonzero elements

This finds all the non-zero points inside the circle with radius 0.4 units, and returns its indices as a column vector.

Next,

u0(ix)=ones(size(ix));

We replace the found non-zero entries of interest by 1 in the zero matrix u0

Next,

tlist=linspace(0,0.1,20);

create an array that will have values of 0 to 0.1 in 20 steps.

Last,

u1=parabolic(u0,tlist,'circleb1',p,e,t,1,0,1,1); % Solve parabolic PDE problem

parabolic is a function that produces the solution to the FEM formulation of the scalar PDE problem. Its input parameters are,

u0 => initial value of the mesh
tlist => time for which solution to be obtained
circleb1 => Initial Boundary condition data
p, e, t => mesh point, edge, triangle data
1, 0, 1, 1 => c, a, f, d
c => viscosity
a => gravity
f => pressure
d => density
u1 => final solution at the end of simulation time.

Now I need to plot this u1.

Joining the Pieces of the Puzzle

12:50 am

So, have been delaying the real implementation so much. Why? I know what all has to be done, how it has to be done, but not sure if it will all work together. So, confirm it all, I have to implement it. But I have been delaying it for so long. Why? I know why.. It’s only because, the work that I have planned is scattered among all the blogs that I have written. Need to put it all together and sit and implement it. Was hoping I could finish all by today, but looks like that’s not happening!! So, today I enlist the things I have to do regarding the model. Here we go.

Step 1 : c₁(p) block

Create a polynomial block with given values of constants, to be multiplied by pressure p which is the input from fetal heart.

Step 2 : c₂(z) block

Create a waveform that approximates the waveform given in the paper. The values could be:

z = [0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80]

c₂(z) = [1 3.5 7.5 10 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11 11.1 11.2 11.3 11.4]

Once that is ready, approximate the input z to the closest value and give out the corresponding c₂(z)

Step 3 : c²(p,z) block

multiply the outputs of step 1 and step 2. This gives a continuous output which is a variable dependent on p.

Step 4 : Location Block

Take the z input and find which arterial segment or bifurcation it belongs to (switch case), and give the distance in that segment as the output.

Step 5 : Initial Cross Section block

Find its initial cross section according to formula for that segment or bifurcation with the output from step 4.

Step 6 : Final Cross Section Block

Multiply the output from step 3 to the blood density constant continuously, find its reciprocal, add it to 1, and then multiply the whole to the output of Step 5. This gives a continuously varying cross section variable.

Step 7 : Stokes – Navier block

Write the S-function which calculates the solution to the stokes navier equation using the Thermal regulator: PDE Toolbox - Simulink application and the following:

Step 8 : Export to the workspace and store, and plot against time

ONLY THE PROBLEM IS PLOTTING AGAINST TIME!!! DISCUSS WITH SIR!!

Confirming the 4th Change in Direction and Some Relief!!

11:05 pm

So, after trying very hard with the method of characteristics used in the original paper, I realized that, the mathematical part is too difficult for a child like me. So, will not be solving the stokes Navier equation using method of characteristics but will be doing it using pdetool in Simulink. This was the same idea that I had posted before. I realize that this means it is a waste of at least a week, but, I’ll call it use,.. since I know, exactly WHY I’m not using the method of characteristics. Have to enlist the reasons.

So, today, before going to sleep, the agenda is, use Simulink to create a module for the Generates the solution to the Stokes-Navier equations. I think the designing part will be easy, but what we need to verify is whether the output we are getting is right or wrong. To do that, we need to connect the entire system together, or give a temporary input that represents the practical input to the module.

1. Think how to calculate mean pressure at any point in the system which I suppose becomes the input to the Stokes-Navier module.

2. Think up if and how this module and the rest of the modules can work together.

3. Form the Stokes-Navier module, using the document mentioned before.

This document will be edited after about half an hour.

So, here I am, trying to design the Stokes – Navier module. Firstly consider this equation given in the document:

Now consider the Stokes-Navier Equation given in one of the ppts that I’m referring to which helps me understand the meaning of the equation:

Consider this Stokes-Navier equation which is given as the prototype in the mathwork manual with good explaination:

Now consider this Stokes-Navier Equation which is given as the prototype in the pdetool with notations which are somewhat different:

Now I have to create an equivalency to match the coefficients. I’m guessing that it’s as follows:

rho = d = density of blood

mu = c = viscosity of blood

The controversy is about a and f. Going by the signs of the terms, the pde toolbox help document and the ppt that I’m referring to, I guess,

p = a = pressure and

f = f = other forces

Now, hoping that these assumptions are correct, I’ve got to substitute the values of these constants. These values are, unfortunately not given in the paper, so I will have to borrow from google, possibly, wikipedia.

Taking, d = 1060 kg/m³, c = (from wikipedia), a will be the time varying mean arterial input pressure, f I think we’ll consider as 1(not sure).

Just figured out how to apply the elasticity to the generated Cross Section of the vessel.

With the Cross section and  space coordinate equations, we can find out a value S of the cross section. Now, as per the formula given in the base paper, we have,

This equation defines elasticity as the steady cross section value, divided by the product of blood density and change in the cross section due to pressure. Now consider the following:

Here, as we know, S becomes the input, we have to have a module for calculation of c².

From here I realise that I have solved the problem of Elasticity application to vessels.

Then, there was this problem about finding the mean pressure.. well just realized how much a of a fool I am to think that pressure in an artery can come from two directions!! Haha.. So, this problem will come only at the placenta. But there too, we have a division as maternal side of the placenta, and the fetal side of the placenta, so, nowhere we will have to encounter a scenario where we will have to calculate the ‘mean’ pressure! HA! Such a fool! Shi.

So now that I have figured out almost the entire deal I am at peace. Aaah! This feels great! There’s a time for everything. Jab jab jo jo hona hai, tab tab so, so, hota hai. Happy!! So, Again! I’ll do all of this tomorrow!!

Now what remains is the thinking up that I have to do regarding the bifurcations.. There is nothing about the bifurcations that is mentioned in the paper! Only the geometrical dimentions are given! What the hell is that all about? I guess, it should be treated in the same way as the arterial segments, but the only difference is that, there are two different values of cross sections at the second end. yoyo! Very happy and relieved! Everything has its own time!

ATB!

Serious Glitches..! Worried Now..

11:36 pm

Have been trying to avoid it, but that's going to be dangerous, gotto face it, and find a way out! Not run away! So here we go!

After reading the "Model of the Main Arterial Tree" part from the paper very carefully I have realised the following things:

The model is based on, firstly, the Mass Conservation Equation given by,

where,

 v is the mean blood flow velocity

S is the arterial cross-section

z is the distance from the aortic valve along the arterial tree

psi is a function describin outflow from the small arteries which are not geometrically depicted

They have mentioned in the paper that they’ve solved this equation using the method of characteristics. Studied it with the help of this video and tried to solve my equation in the following way: (hope it is all understandable)

As you can see, I encounter many problems while trying to solve the mass conservation equation.

Firstly, I do not know if v is a variable which is dependent on the value of z or not. Thinking about it, we realize that the velocity does actually change depending on its distance away from the heart. But that happens only due to the change in pressure with distance which we rae going to accommodate for. If we include this dependence, the equation becomes difficult to be solved using the method of characteristics. Saying this, coz I tried doing that as shown below:

Thus, I decided to assume it independent of z, and continued as shown in the first image.

Secondly, knowing that v is the mean velocity, we try to understand the meaning of the word mean velocity. Considering the fact that in the paper, they have mentioned that the objective of using this method is to be able to compute blood velocity and pressure at any point in the arterial tree as the contribution of the incident and reflected velocity and pressure. Which implies that, by mean velocity, they mean, the resultant of velocity and pressure from maternal side and the fetal side. This makes us think that the v mentioned, needs to be the input to the simulink block that will solve the mass conservation equation, and before entering this block, there needs to exist a simulink block that will calculate this mean velocity. At whichever point in the arterial tree we  need to check the velocity or pressure, we need to do the calculation for that part.

This conclusion gives rise to a controversy. Is this input velocity the same that will be shown as the output velocity? Cannot accurately answer this question right now, coz we haven’t come to that point in the process still, but all that I can guess right now, is that, this input velocity will be acted upon by the mass conservation equation, then the stokes navier equation etc. and get modified. Coming to think of it now, I feel that the input to the system is this mean velocity, but the output might be the pressure wave. This might be needed to be converted into velocity somehow, since none of the equations in this part give velocity as the output.<== Problem

Going into the depth of the above mentioned problem, now we need to

1. Figure out how to calculate the mean velocity
2. Find the expression that describes the relation between velocity and pressure of an incompressible fluid.

Once the mean velocity calculation problem is solved, we need to know how to incorporate this into the mass conservation equation, and design a simulink block that implements the mass conservation equation.

So now, let us talk about this solution of the mass conservation equation, In the equation, they have mentioned psi, whose value hasn’t been mentioned anywhere in the paper. Need to ask what value to input as PSI. Assuming currently that it is zero, I move ahead. As shown in the first image, this gives the solution as

S= f(z_o)=f(z-vt)

This equation says that the cross sectional area for the curve represented by z-vt = z_o is given by the function which is the initial condition at t=0 for that curve.

For all arterial segment we have been given the length of the segment and the cross sectional area at its 2 ends. The table is as follows:

Now as the solution the cross section needs to be a function of z-vt. What function? As the tutorial video suggests, this function is the initial condition given. What is the initial condition given to us? These are just numbers. Thus we need to form relationships, preferably, linear from the given values of cross sections and length. To do this, I used curve fitting tool from MATLAB. Here’s the video of how I did it. Sorry for the delay between video and audio, and the overall bad quality.

Okay, the video refuses to get uploaded, sorry! Here are the equations that I got:

For Arterial segment 2,   

For Arterial segment 3,   

For Arterial segment 4,  

For all the above expressions, y is the respective Cross section and x is the respective position coordinate.

For arterial segment 1 and arterial segment 5, the cross section remains constant. That is it does not change for any value of length. Thus the initial value remains the same for all space coordinates.

Now that we have a relationship that describes the initial condition and relates cross section with the z coordinates, we need to figure out how to incorporate the value of z into the initial condition and how to incorporate the initial condition into the solution.

For each of the above x, we substitute z- vt where z is the distance from the aortic valve along the arterial tree to the beginning / end (?) of the arterial segment, v is the computed mean velocity and t is the simulation time. So now, we have got the logic for doing a part of the Arterial segments modeling.

After this, is the block that computes the arterial vascular elasticity and incorporates that into the model. This part is given as follows:

So, c₁(p) is a polynomial block, and c₂(z) is a waveform to be approximated. these can be implemented on SIMULINK.

All this to be implemented on MATLAB and have to think about the mean velocity block, tomorrow!

ATB!

Hello World!

04:43 pm

Trying to work on the project, but my cousin sister keeps disturbing me.. Making aweful noises in my ear.. Cheh.. But hey!! It’s fun!!

Trying out Windows Writer for the first time.. Let’s see how this goes!

Yo!! ATB!

Edit1: My sis angry! So, … Credit of the Writer Discovery to her!!! Maskaa!!

Modeling the Feto-Maternal Circulation Arterial Tree

07:52 pm

Last 2 Sundays have been spent in preparing the rough sketch of discrete parts of the model. To begin with, I divided the model into:

1. Fetal Heart (Pressure Source 1)

2. Maternal Uterine Arteries (Pressure Source 2)

3. Peripheral areas:
(PA.1) Brain
(PA.2) Kidneys
(PA.3) Peripherals at Primitive Iliac Artery Bifurcation
(PA.4) Peripherals at External Iliac Artery Bifurcation
(PA.5) Peripherals at Internal Iliac Artery Bifurcation

4. Arterial Segments:
(Seg.1) Ascending Aorta (Ductus Arteriosus)
(Seg.2) Thoracic Aorta
(Seg.3) Abdominal Aorta
(Seg.4) Primitive Iliac Artery
(Seg.5) Internal Iliac Artery
(Seg.6) Umbilical Arteries
(Seg.7) Uterine Arteries

5. Bifurcations:
(Bif.1) Cerebral Arteries
(Bif.2). Renal Arteries
(Bif.3) Primitive Iliac Arteries
(Bif.4) External Iliac Arteries
(Bif.5) Internal Iliac Arteries

6. Placenta on the fetal side of the feto-maternal circulation

7. Placenta on the maternal side of the feto-maternal circulation

1. Modeling of the fetal heart:

Fetal Heart parameters, ventricular pressure and ventricular volume follow the relation,

p(t) = E(t)[V(t) - V0]

where,
p(t) => ventricular pressure in mmHg
V(t) => ventricular volume in ml
V0 => Reference Volume
E(t) => Ventricle Elastance

In this relation, p(t) and V(t) are unknowns. Thus, I had to create a pressure wave that approximates the one given in the paper. I used the following values:

time axis => [0, 0.049, 0.054, 0.059, 0.07, 0.08, 0.1, 0.149, 0.173, 0.198, 0.208, 0.217, 0.225, 0.251, 0.4]
pressure axis => [0, 0, 1.09, 1.49, 1.79, 1.93, 1.99, 1.99164, 2, 1.96, 1.9, 1.82, 1.65, 0, 0]

To find E(t), there exists an equation

E(t) = Emax En(t)

where En(t) = 5.412t^6 - 20.066t^5 + 25.542t^4 - 13.71t^3 + 2.714t^2 + 1.08t + 0.029

Here t is taken as the simulation time. Thus, to create the block, now we need to have the values of Emax and V0.

After the analysis of fetal lamb and dog heart and approximately equating it to human fetus, we obtain, Emax = 6 mmHg per ml and V0 = -8ml.

Error: The output produced shows pulsetile increase in pressure, but it does not reduce to the baseline!

2. Modeling of the maternal uterine arteries:

I approximated the above shown waveform which was the only information provided about the uterine pressure. The values that I used are:

Pressure values => [78, 77, 78, 80, 90, 100, 110, 118, 122, 130, 131, 130, 122, 118, 112, 108, 103, 100, 90, 87, 84, 82, 84, 84, 80, 79, 78]
Time Values => [0, 0.05, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22, 0.24, 0.25, 0.26, 0.28, 0.3, 0.32, 0.34, 0.36, ,0.38, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8]

3. Peripheral Areas

Q(t) = {[p(t) - pc] / R} + {C dp(t)/dt}

where
Q(t) => Blood flow
R => Resistance of PA
C => Compliance of PA
pc => Critical Pressure of PA

With the above mentioned values from the paper and the equation mentioned, I created the PA basic subsystem in Simulink as shown:

and varied the pc, R, and C values for every PA.

Arterial Segments and Bifurcations : Main arterial tree

1. In the paper that we are referring to, they have specified the geometrical parameters on the fetal side of feto-maternal arterial tree for each of the segment and bifurcation for fetus at term. BUT, In reference to the blood flow in these segments they have NOT given any values for the parameters like

• the density of blood 
• viscosity 

which are necessary in the formation of the Navier-Stokes equation. This equation, together with mass and momentum conservation equationsform the constitutive equations of the model.

Suggested Solution: We could use values from any other reference.

2. This highly mathematical part of the model in the paper, which requires the solution to the Navier - Stokes Equation as solved using Method of Characteristics, has not been provided.

Suggested Solution: We could use the pdetool in MATLAB to solve the Navier- Stokes Equation for incompressible fluid. For this, I found this attached document on the internet. Asked my External Guide to go through it and see if we could use it, He happily showed the green signal. Now gotta move onto the next question.
Since we see that in the above mentioned paper, twhy use the PDETool, to solve the fluid flow equations, and I am using Simulink, I need to somehow interface these two. For this I just found this paper very useful.It gives an example in which they have attempted an interface between, Simulink and PDETool. Will be studying and implementing the model from these sources!

ATB!

Mathematical modeling, glitches expected

05:35pm

Yes, I am back after a long,.. long tiime! Yes, been bus, No, not studies.., watching movies.. :) Not really though. Internal tests were there. Had to take a break from project work. So here we begin again. Discussed various models with the External Guide, and finally came up with this paper as the most ideal one.

Here's the jist of what was discussed on the last callcon with my External Guide:

The selected model is an ideal model because:

1. It simple.
2. Model parameter values have been provided.
3. It not only talks about feto-placental circulation, but gives circulation in other important arteries like cerebral artery.

Short Timeline:

• By 2 weeks, that is, till Diwali, the simulation of this model on Simulink should be ready for review by the External Guide.
• If the model gives results that match the ones given in the paper, for a normal lady, then pathological cases could be taken up.
• This could be then presented to the ethical committee, and approved for collection of data for abnormal patients.
• Parallely, could think about writing a paper. 
• How we are going to present this topic?
• Which conference to select? (Preference: International, IEEE conference)

PS: The experimental setup that the HOD was talking about should not be necessary anymore.