Analyses and Simulation of Medicine Carrying Blood Flow in MATLAB

1. Analyses and Simulation of Medicine Carrying Blood Flow in MATLAB MTH 499 Jill Mercik and Ryan Banci

2. Our Research Project • Simulating through MATLAB how blood flows through the human circulatory system. • 1D simulations • Advisor: Prof. Yanlai Chen

3. Why are we doing this? • Simulations of the natural flow of blood can give us information on how different medicines are distributed through the body. • This information can help researchers to design specific treatments that can improve upon patient care.

4. Overview • Circulatory System • Mathematical Techniques • Finite Difference Method • Navier-Stokes Equation • Numerical Scheme • MATLAB • 1D Simulations • Conclusion • Future Work

5. Circulatory System Made up of: • Arteries • Arterioles • Capillaries • Veins • Heart Blood is made up of: • Cells • plasma

6. Mathematical Techniques Finite Difference Method • Used to approximate the solutions to differential equations. • Algebraic in form and solutions are related to grid points.

7. Finite Difference Method Example Consider the ODE: The Euler method uses the finite difference quotient: In order to approximate the differential equation by first substituting and applying some algebra to get: Solving this last equation can give us an approximate solution to the differential equation.

8. Finite Difference Method We will use implicit finite difference. Using the given equations we will create a system of equations that we will put in matrix form in order to be solved in MATLAB. The solutions of this matrix of equations will give us the values for A(cross sectional area) and q(flow rate) that we need to simulate the branch for time domain n+1. We can then apply this solution to solve the next sequential time domain then repeat until we have reached our requested simulation time. How it is used in our research project

9. Mathematical Techniques What is it?? Nonlinear PDE’s the describe the motion of fluid substances. Comes from applying Newton’s 2nd law to fluid motion. F = ma The Navier-Stokes Equations are used in a lot of practical problems such a modeling the weather, ocean current, and water flow in a pipe. Navier Stokes Equations

10. Navier-Stokes Equations General form of the Navier-Stokes equations that are used to model blood flow are expressed as: Equation (1) is a form of the momentum equation while equation (2) is a volume continuity equation that was simplified from the mass continuity equation.

11. Navier-Stokes Equations The previous slide’s equations can be written in component form as:

12. Numerical Scheme In order to simulate a branch of the blood stream, we need to find the flow rate (q) and cross-sectional area (A) for a section of points along the branch and repeat this for each time step. The more points (k) we use, the more accurate the calculation results will be.

13. Numerical Scheme We need a few equations in order to make a 1D simulation. The equations for flow and pressure can be expressed through the conservation of mass and momentum as:

14. Numerical Scheme From those equations, we can derive the equations we need to find a system of equations that we will put into matrix form to be solved in MATLAB. We found diagonal patterns in the matrices that we could code in a MATLAB script file to generate automatically given a user input for the sample size k.

15. MATLAB

16. MATLAB %Variables; k=input('The sample size is: '); dt=input('The change of time (delta t) is: '); dx=input('The change of position (delta x) is: '); K =(dt/dx); E =(dt/(dx^2)); q = 1; %(68*10^-3); A = 1; %(2*10^-11); p = 1; h0 = (6000*10^-9); R0 = (2209*10^-9); A0 = (1.5333*10^-11); v = 1.1; Estat = (7*10^-6); a = (8/3)*(q/A); b = ((-4/3)*(q/A)^2)+((A/p)*((Estat*h0)/(R0*A0))); d = v;

17. MATLAB %Inlet conditions (B4, B5, B6) g1=[1 zeros(1,k-1) (K/2) zeros(1,(k-2))]; g2=[0 1 zeros(1,k-2) 0 (K/4) zeros(1,(k-3))]; g3=[(-b*K/4) 0 (b*K/4) zeros(1,k-3) (1+E*d) (a*K/4)-(E*d/2) zeros(1,(k-3))]; %Outlet conditions (B12, B13) g4=[zeros(1,k) 1 zeros(1,k-5) (-K/2) 0 (K/2)]; g5=[zeros(1,2*k-2) 1];

18. MATLAB %Middle Matrix Begins (B2, B3) x1 = ones(1,k-3); x2 = diag(x1,0); y1 = (K/4)*ones(1,k-3-1); y2 = diag(y1,1); y3 = (-K/4)*ones(1,k-3-1); y4 = diag(y3,-1); yout = y2 + y4; j1 = (b*K/4)*ones(1,k-3-1); j2 = diag(j1,1); j3 = (-b*K/4)*ones(1,k-3-1); j4 = diag(j3,-1); jout= j2 + j4; z1 = ((a*K/4)-(E*d/2))*ones(1,k-3); z2 = diag(z1,0); z3 = (1+E*d)*ones(1,k-3-1); z4 = diag(z3,-1); z5 = ((-a*K/4)-(E*d/2))*ones(1,k-3-2); z6 = diag(z5,-2); zout = z2 + z4 + z6; b1 = zeros(2*k-6,2); %buffer zeros b2 = zeros(2*k-6,3); %Middle Matrix Ends

19. MATLAB

20. 1D Simulations • After solving the matrices with the given boundaries, we are able to simulate the branch by applying the values of A (Cross Sectional Area) and q (Flow rate) to each part of the branch for one time domain. Then using the data from that time domain to solve the next sequential domain.

21. Conclusion This project ended up being a lot harder than we originally thought. We had hoped to create some 2D simulations, but were only able to create 1D simulations. Modeling blood flow is a lot harder than modeling something like water in a pipe because of the different pulse flows and internal/external forces. Future Work We need to apply the program to a large sample size of 100 or greater in order to solve for the n+1 time domain. Then we would analyze the plots for flow rate vs. time and for flow vs. area. Also, we will use the solutions for the solved n+1 time domain to solve the next sequential time domains. This would complete our 1D simulations. Then, the next logical step would be to create 2D simulations, and ultimately 3D simulations.