Intracranial Dynamicsbioe 310 - Biological Systems Analysis

Intracranial DynamicsBioE 310 - Biological Systems Analysis

Intracranial DynamicsandHarmonic Analysis

Richard Hickey

As part of the coursework for BioE310

Report produced under supervision of Professor AndreasLinninger

Dept. of Bioengineering, University of Illinois at Chicago

Abstract

A quantitative model of blood flow in intracranial vasculature, in tandem with cerebrospinal fluid (CSF) flow in the intrathecal space, is created to better understand and predict how pressure changes in one system will affect the other.A compartmental approach is used to analyzeflow and pressure changes between major areas of the intracranial space. Finally a harmonic analysis provides insight into the role of arterial pulsation in whole-brain fluid dynamics at ultimate periodic response.We find that the amplitude and phase angle of the pulse differ between each compartment, illustrating how the various compartments all react uniquely to each heartbeat.

12/12/141

Intracranial DynamicsBioE 310 - Biological Systems Analysis

1. Introduction

The dynamic nature of fluid flow throughout the cranium is of significant interest to many disciplines of medicine. The minute pressure changes within the intracranial space that occur during each cardiac cycle (heartbeat) have profound effects on the efficacy of intrathecal drug delivery[2]. The relationship between arterial, capillary, and venous pressures within the body are well studied. However there are many compartments within the cranium that are difficult to assess directly. The activity within these areasis of paramount importance when studying pathological changes in intracranial pressure (ICP) such as in hydrocephalus or traumatic brain injury[7]. Despite a demand for detailed information of this domain, the invasive measures currently required to obtain such parameters are contraindicated in all but the most extreme cases[5,7].

Given that the system can be measured at the inlet (carotid arteries) and outlet (jugular veins), it is reasonable to attempt to construct a model of what activity takes place within the intracranial system. The arterial pulse can be treated as a waveform and as a sort of “forcing function” that propagates activity all through the brain and surrounding compartments. The study of how this one waveform interacts to stimulate other similar, but not identical, waveforms is called harmonic analysis. It is through an understanding of these harmonics that we can model the activity in areas not accessible for direct instrumentation.

2. Methods

This modelworks with the 6 major compartments of intracranial space: Cerebral arteries, capillaries, veins and venous sinus, as well as brain and CSF[5,6]. Limited numericdata regarding intracranial flowis available. Some of the more comprehensive work, done by Sorek et al,provides constants for resistance and compliance of the various intracranial compartments[5]. These values were used in a matrix of coefficients that describe a mass conservation balance equation. Boundary conditions were obtained using the Monro-Kellie hypothesis of the fixed-volume cranial space and standard assumptions of mean arterial pressure[1,4]. Pressures and flows of the steady state system were solved using MATLAB matrix operations.

The harmonic analysis was done with single-variable forcing equation (1), and afrequency of 80bpm was chosento roughly approximate an arterial waveform[8].

/ (1)

Finally, pressures from the earlier computation were introduced into the harmonic system to create a more comprehensive simulation of intracranial dynamics. A 5-variable, single-frequency, first order system was input into MATLAB in the form of a 5x5 matrix Z.

/ (2)
/ (3)

When f0 is a column vector consisting of values for I0 (from the forcing function) and 0, a solution X exists with complex entries that describe the harmonic waves.

3. Results

A MATLAB script was used to generate a map of the intracranial space. This flow diagram indicates which compartments communicate with one another. Solutions for mean pressure and flow from the steady state equation are plotted on the diagram (Fig. 1). The inflow/outflow calculated are consistent with literature values at approximately 750ml/min[3,6].

Figure 1: (Top) Conceptual flow network showing communication between compartments
(Bottom) Equivalent MATLAB generated plot, with flows (Q) and pressures (P) labeled

3.1 Pressure Model

Using the same set of equations plus a periodic function for arterial waveform, a simple graph of pressures/time can be visualized (Fig. 2). This is useful for observing trends in pressure and is consistent with the work of Zhou, Xenos et al[8].Each of the 6 compartments can be seen oscillating at different levels of pressure. However, this model (like [8]) isincomplete because it does not account for any damping of the waveform.As it travels through a significant amount of vasculature and tissue around the brain, and as the mean pressure drops, it would be expected that the amplitude of the lower pressure waves would be lower than the initial arterial one. Indeed that is the case, which is something we can observe using the more thorough mathematical modeling that is harmonic analysis.

Figure 2: The arterial waveform reproduced at steady state pressures for all compartments

3.2 Harmonic Analysis

Harmonic analysis was applied to the sinusoidal function representing arterial flow using 5x5 matrixZ as in equation (2). Because of the sine wave representing arterial flow (1), only the imaginary part of X is preserved. The results of one complete oscillation, omitting the baseline pressures, are shown in (Fig. 3) and demonstrate the phase shift between compartments. These harmonic waves are all shown at the ultimate periodic response and have the same period as the arterial pulse. They can also be visualized individually (Fig. 4).

Figure 3: Over a single cycle of arterial flow (black), each compartment can be seen with a unique amplitude and phase shift.

These diagrams are made using the assumption that the heartbeat is a continuous signal (we hope!) and so does not include the introductory period in which the responding waves are still adapting to the forcing (arterial) function. Rather they are shown at their ultimate periodic response, as time t becomes arbitrarily large.

Although the waves still appear to vary greatly, they all share the same frequency as the initial (arterial) wave. Where the waves differ is in two values: their amplitude and their phase shift. The baseline for each measurement comes from the arterial wave, as that is the forcing function. All other waves are shown as percentage of arterial pressure and delay vs. arterial frequency. There is a predictable stepwise drop in pressure as the fluid moves from artery to capillary to vein. Interestingly, the brain compartment experiences a markedly different phase from the other compartments. A working theory is that because it is being fed by more than one inlet (see flow diagram inFig. 1) its waveform has a pattern unique from the others.

Figure 4: (Top)Harmonic waves of individual compartments influenced by the arterial wave at 80bpm(Bottom) Amplitude and phase shifts relative to the arterial pulse.

Figure 5: Pressure values are combined with the harmonic waves to create a comprehensive plot of intracranial activity

With this better understanding of the waveforms as they exist in a real system with resistance, we can revisit the earlier plot of pressure waveforms vs. time. This new analysis represents our most complete model, which incorporates mean pressures as well as harmonics (Fig. 5).

4. Validation

The results seen in this study accurately predict those seen in a selection of other works [6,8]. The shift in amplitude observed inFig. 2 reflect reasonable physiologic differences between intracranial compartments.

Figure 6: Reproduced from Zhou, Zenos, et al [8]. Generally similar pressures were observed in both systems. However this image, like Fig. 2, lacks a complete harmonic adjustment of wave amplitude.

5. Discussion

These preliminary results suggest a significantly pulsatile nature of pressure within the brain tissue itself, which if accurate has interesting implications warranting further study. Such a finding is consistent with some previous literature including the work of Wagshul et al and leads to a more thorough understanding of how changes in blood pressure affect the pulsation of CSF[7].

A shortcoming of this model is the simplistic waveform chosen to represent the cerebral artery pulse. A single sine wave was used because its harmonic influence on a multivariable system is well understood. In future studies it would be of interest to perform a Fourier transform of a measured physiological signal in order to obtain a more accurate representation of the pulse waveform. A similar model could then be created using this multiple-frequency forcing function. This would make for a more complete model but no significant change in results would be anticipated.

One major assumption made in this study is that flow into the intracranial space is always equal to flow out. Aspects of this Monro-Kellie theory have been upheld for over 2 centuries [4]. Principally that due to the rigid bony structure of the skull there must be a fixed volume inside of it at all times. The balance of this volume between brain, CSF, blood, and other matter may vary but the total volume must not. A caveat that has not been properly examined by subsequent “compartmental” type approaches is that not all of the compartments reside fully within the intracranial space. The major culprit, of course, is the CSF. With each cardiac cycle a significant amount of CSF (2ml) pulsates down into the spinal cord. The spine does not follow the same set of rules regarding compressibility – namely, the vertebrae do not provide the same level of confinement as does the skull. A modified flow network that ventures outside of the intracranial space and into the spine is outlined below (Fig. 7). In future studies it would be beneficial to examine what sort of impact the detour of CSF down the spinal column has on the rest of the intracranial network.

Figure 7. Is there another compartment not adequately represented? The spinal cord is not within the intracranial space but communicates with it through the CSF.

6. Conclusion

It is possible to create a compartmental model of the intracranial system to mimic what occurs inside that difficult-to-access part of the human body. The combination in this model of pressures and harmonics appears to represent a novel and more complete simulation of intracranial flow. Much previous work in the field had not accounted for signal damping in the resistive tissue in the intracranial system[8]. Based on the results from the combined pressure & harmonics model (Fig. 5) we can see that the waveform within intracranial structures varies greatly.

This model demonstrates the power of harmonic analysis to describe the inner workings of what is otherwise an inaccessible system. These findings may be used to better design systems for intrathecal drug delivery or to better understand pathological changes that occur with rising ICP.

7. Acknowledgments

Report produced under the supervision of Bioengineering Professor Andreas Linninger, with the generous aid of Chih-Yang Hsu and Sebastian Pernal.

Intellectual Property: Biological and physiological data and some modeling procedures provided to from Dr.Linninger’s lab are subject to IRB review procedures and intellectual property procedures.Therefore, the use of these data and procedures are limited to the coursework only. Publications need to be approved and require joint authorship with staff of Dr.Linninger’s lab.

8. References

1.Dixon, W., and W. Halliburton. The cerebro-spinal fluid. II. Cerebro-spinal pressure. J. Physiol. , 1914.at <

2.Hsu, Y., H. D. M. Hettiarachchi, D. C. Zhu, and A. Linninger. The frequency and magnitude of cerebrospinal fluid pulsations influence intrathecal drug distribution: key factors for interpatient variability. Anesth. Analg. 115:386–94, 2012.

3.Lakin, W., and S. Stevens. A whole-body mathematical model for intracranial pressure dynamics. J. Math. Biol. 383:347–383, 2003.

4.Mokri, B. The Monro-Kellie hypothesis: Applications in CSF volume depletion. Neurology 56:1746–1748, 2001.

5.Sorek, S., J. Bear, and Z. Karni. A Non-Steady Compartmental Flow Model of the Cerebrovascular System. Biomechanics 21:695–704, 1988.

6.Stevens, S. A. Mean Pressures and Flows in the Human Intracranial System, Determined by Mathematical Simulations of a Steady-State Infusion Test. Neurol. Res. 809–814, 2000.

7.Wagshul, M. E., P. K. Eide, and J. R. Madsen. The pulsating brain: A review of experimental and clinical studies of intracranial pulsatility. Fluids Barriers CNS 8:5, 2011.

8.Zhou, X., M. Xenos, and S. Kontapalli. Response on Harmonic Excitation Analysis. Laboratory for Product and Process Design, University of Illinois at Chicago: 2005.

12/12/141

Intracranial DynamicsBioE 310 - Biological Systems Analysis

Appendix A: IntracranialPressure.m

12/12/141

Intracranial DynamicsBioE 310 - Biological Systems Analysis

clear all; close all; clc;

%% Createvectors

% Positions to graph flow network

ptCoordMx = [5 10; 5 7; 2.4 6.3; 7.5 6; 5 4; 5 3; 5 2;];

%Constitutive relationships

pointMx = [-1 0 0 0;

1 -2 -5 -3;

2 -4 -7 0;

3 4 -6 0;

5 6 -8 0;

7 8 -9 0;

9 0 0 0];

faceMx = [1 2; 2 3; 2 4; 3 4; 2 5; 4 5; 3 6; 5 6; 6 7];

%% Given values

Rac = 0.106667; %mmHg s / ml

Rcv = 0.020;

Rbv = 0.6;

Rvs = 0.0027;

Rfs = 2.80;

Rcb = 1.4;

Rfb = 10.00;

Rcf = 1;

alpha = [Rac Rcf Rcb Rfb Rcv Rbv Rfs Rvs 0];

[Plength,Pwidth]=size(pointMx);

[Flength, Fwidth] = size(faceMx);

Findex = 1:Flength;

Pindex = (Flength+1):(Flength+Plength);

Matindex = [Findex Pindex];

Matlength = length(Matindex);

b = zeros(Matlength,1);

k = 1;

Pin=1;

Pout=7;

b(Pin)=102;

b(Pout)=5;

t=0:.01:5;

I0 = 6; % amplitude

HR = 80; % in BPM

w0 = 2*pi*(HR/60);

Pa = 102 + I0*sin(w0*t);% Arterial waveform

Ps = 5 + I0*sin(w0*t);% Outflow waveform

for T =1:length(t)% BEGIN TIME LOOP

k=1;

C = zeros(Matlength,3); %C is column vectors Row, Col, Value

b(Pin) = Pa(T);

b(Pout) = Ps(T);

%% Conservation Eqs: F1 - F2 = 0

for i=1:(Plength)

%Exception: first

if(i==Pin)

C(k,1) = i;

C(k,2) = Pindex(i);

C(k,3) = 1;

k=k+1;

%Exception: last

elseif(i==Pout)

C(k,1) = i;

C(k,2) = Pindex(i);

C(k,3) = 1;

k=k+1;

else%NOT first or last

for j=1:Pwidth

if pointMx(i,j) < 0

C(k,1) = i;

C(k,2) = abs(pointMx(i,j));

C(k,3) = -1;

k=k+1;

end

if pointMx(i,j) > 0

C(k,1) = i;

C(k,2) = abs(pointMx(i,j));

C(k,3) = 1;

k=k+1;

end

end

end

end

%% Constitutive Eqs: a1F1 - P1 + P2 = 0

for i=1:Flength

for j=1:Fwidth

if j==1

C(k,1) = i+Plength;

C(k,2) = abs(faceMx(i,j)+Flength);

C(k,3) = -1;

k=k+1;

end

if j==2

C(k,1) = i+Plength;

C(k,2) = abs(faceMx(i,j))+Flength;

C(k,3) = 1;

k=k+1;

end

end

C(k,1) = i+Plength;

C(k,2) = i;

C(k,3) = alpha(i);

k=k+1;

end

% Populate SPARSE Matrix

B = sparse(C(:,1),C(:,2),C(:,3));

% Determinant and solution

d = det(B);

x = B\b;

Yplot(T,:)=x(Pindex);

end %END TIME LOOP

%% Plot Pressures

plot(t,Yplot(:,1:end-1));

axis([0 5 0 112])

legend('Pa','Pc','Pf','Pb','Pv','Ps','Location','East');

%% Create Network Diagram

figure

hold on

for i=1:length(ptCoordMx)

scatter(ptCoordMx(i,1),ptCoordMx(i,2));

pointS = sprintf(' P%i=%.1f',i,x(Pindex(i)));

text(ptCoordMx(i,1),ptCoordMx(i,2),pointS);

end

for j=1:Flength

plot([ptCoordMx(faceMx(j,1),1) ptCoordMx(faceMx(j,2),1)],...

[ptCoordMx(faceMx(j,1),2) ptCoordMx(faceMx(j,2),2)]);

xspot = (ptCoordMx(faceMx(j,1),1)+ptCoordMx(faceMx(j,2),1))/2;

yspot = (ptCoordMx(faceMx(j,1),2)+ptCoordMx(faceMx(j,2),2))/2;

pointS = sprintf(' Q%i=%.1f',j,x(Findex(j)));

text(xspot,yspot,pointS);

end

12/12/141

Intracranial DynamicsBioE 310 - Biological Systems Analysis

Appendix B: Harmonics.m

12/12/141

Intracranial DynamicsBioE 310 - Biological Systems Analysis

clear all; close all; clc;

% P0 values from IntracranialPressures.m

P0 =[102;%a

21.8122607460216;%c

16.9962366130515;%f

11.6796974728199;%b

7.01817800731953;%v

5];%s

% Constants, from SOREK et al

Cab = 0.0155;

Ccf = 0.0364;

Cbv = 0.3180;

Cfs = 0.0334;

Cfb = 0.1830;

Rac = 0.106667; %mmHg s / ml

Rcv = 0.020;

Rbv = 0.6;

Rvs = 0.0027;

Rfs = 2.80;

Rcb = 1.4;

Rfb = 10.00;

I0 = 6; % amplitude

HR = 80; % rate in BPM

w0 = 2*pi*(HR/60); % angular frequency (rad)

% Time to plot (s)

% t=0:.01:(60/HR); % Single cycle

t = 0:.01:5;% 5 second plot

% (a) c f b v s

Z = [ ...

1i*w0*Ccf+1/Rac 0 0 0 0;

0 1i*w0*Cfb+1/Rcv 0 0 0;

0 0 1i*w0*(Cbv-Cab)+1/Rbv 0 0;

0 0 0 1i*w0*(-Cbv)+1/Rvs 0;

0 0 0 0 1i*w0*(-Cfs)+1/Rvs];

A = [real(Z) -imag(Z);

imag(Z) real(Z)];

y = [I0 I0 I0 I0 I0 0 0 0 0 0]';

X = A\y;

Xc = [ ...

X(1) + 1i*X(6);

X(2) + 1i*X(7);

X(3) + 1i*X(8);

X(4) + 1i*X(9);

X(5) + 1i*X(10)];

for n=1:length(t)

V(n,:) = imag(Xc.*exp(1i*w0*t(n))); % take only imaginary part

end

Y = I0 * sin(w0*t); %original signal

plot(t,Y+P0(1),'k',t,V(:,1)+P0(2),':r',t,V(:,2)+P0(3),'-.m',t,V(:,3)+P0(4),'--b',t,V(:,4)+P0(5),'-g',t,V(:,5)+P0(6),'--y')

legend('1: Arteries','2: Capillaries','3: CSF','4: Brain','5: Veins','6: Venous Sinus','Location','East')

xlabel('Time (s)'); ylabel('Pressure (mmHg)');

title('Allcompartments, with initial values');

grid on;

figure;

subplot(3,2,1);

plot(t,Y+P0(1),'k');

title('1: Arteries');

subplot(3,2,2);

plot(t,V(:,1)+P0(2),':r');

title('2: Capillaries');

subplot(3,2,3);

plot(t,V(:,2)+P0(3),'-.m');

title('3: CSF');

subplot(3,2,4);

plot(t,V(:,3)+P0(4),'--b');

title('4: Brain');

subplot(3,2,5);

plot(t,V(:,4)+P0(5),'--g');

title('5: Veins');

subplot(3,2,6);

plot(t,V(:,5)+P0(6),'--y');

title('6: Venous Sinus');

12/12/141