# Giant low-density lipoprotein (LDL) accumulation in multi-layer artery wall models¶

In this document the source code of LDL transport simulations. We solve numerically the transport equations:

\begin{eqnarray} \frac{\partial c}{\partial t} +(1-\sigma)\vec u\cdot\nabla c & = & D_{eff} \nabla^2 c - k c \end{eqnarray}

in four layers.

In [1]:
%pylab inline

import numpy as np
from scipy.sparse import dia_matrix
import scipy as sp
import scipy.sparse
import scipy.sparse.linalg

import matplotlib
import matplotlib.pyplot as plt

newparams = { 'savefig.dpi': 100, 'figure.figsize': (12/2., 5/2.) }
plt.rcParams.update(newparams)

params = {'legend.fontsize': 8,
'legend.linewidth': 0.2}

plt.rcParams.update(params)

Populating the interactive namespace from numpy and matplotlib


## WSS dependent parameters¶

The impact of WSS on the transport properties is given by the equations (8-11) developed by Olgac et al.

In [2]:
def phi(WSS = 1.79):
Rcell = 15e-3 # mm
area=.64 # mm^2
SI = 0.38*np.exp(-0.79*WSS) + 0.225*np.exp(-0.043*WSS)
MC = 0.003797* np.exp(14.75*SI)
LC = 0.307 + 0.805 * MC
phi = (LC*np.pi*Rcell**2)/(area)
return( phi)

def Klj(w=14.3e-6,phi=5e-4):
"permability  in m^2"
Rcell = 15e-3 # mm
return ( (w**2/3.)*(4.*w*phi)/Rcell * (1e-6) )

def Kend(w=14.3e-6,phi=5e-4):
"permability  w m^2"
Kend_70mmHg =3.22e-21
Knj = Kend_70mmHg - Klj() # at 70mmHg
return Knj + Klj(w,phi)

def sigma_end(phi=5e-4,w=14.3*1e-6,r_m = 11e-6):
a = r_m/w
Kend_70mmHg =3.22e-21
Knj = Kend_70mmHg - Klj() # at 70mmHg
sigma_lj = 1-(1-3/2.*a**2+0.5*a**3)*(1-1/3.*a**2)
return  1 - ((1-sigma_lj)*Klj(phi=phi))/(Knj+Klj(phi=phi))

def Diffusivity(w=14.3e-6,phi=5e-4,r_m = 11e-6):
"Diffusivity  w um^2/s"
R_cell = 15e-6 # m
a=r_m/w
D_lumen=2.71e-11
return D_lumen*(1-a)*(1.-1.004*a+0.418*a**3-0.16*a**5)*4*w/R_cell*phi*1e-3*1e12


## Parameters¶

In simulation four sets of parameters could be used. Three of them corresond to four layer model. The last one is two layer model.

In [3]:
class LDL_Parameters_Vafai2012(object):
""" S. Chung, K. Vafai, International Journal of Biomechanics 45(2012)"""
names       = [  'endothel' , 'intima', 'IEL'   ,'media' ]
D           = [   5.7e-12     , 5.4e-6  , 3.18e-9 , 5e-8   ]
V           = [   2.3e-2 ]*4
sigma       = [   0.9888    , 0.8272  , 0.9827  , 0.8836 ]
L           = [   2.        , 10.     ,    2.   , 200.   ]
k_react     = [   0.        , 0.      ,    0.   , 3.197e-4 ]
K           = [   3.22e-15  , 2e-10   ,4.392e-13, 2e-12   ]
mu          = [   0.72e-3   , 0.72e-3 , 0.72e-3 , 0.72e-3 ]

def calculate_filration(self,dPmmHg):
mmHg2Pa = 133.3
dP = mmHg2Pa*dPmmHg
Rw = [L_*mu_/K_ for L_,K_,mu_ in zip(self.L,self.K,self.mu)]
self.Vfiltr = dP/sum(Rw)

def __init__(self,WSS=1.79):
""" Change units to mikrometers
Class can be initialized with a value of WSS in Pa
"""
dPmmHg=70
self.phi = phi(WSS=WSS)
self.sigma_end = sigma_end(phi=self.phi,w=14.3*1e-6,r_m = 11e-6)
self.Kend = Kend(w=14.3e-6,phi=self.phi)
self.K[0] = self.Kend*1e6
self.calculate_filration(dPmmHg=dPmmHg)

self.D = [D_*1e6 for D_ in self.D]
self.V = [ self.Vfiltr*1e6]*4
self.sigma[0] = self.sigma_end
self.D[0]=Diffusivity(w=14.3e-6,phi=self.phi,r_m = 11e-6)
def get_params(self):
return (self.D,self.V,self.sigma,self.L,self.k_react)

class LDL_Parameters_Vafai2006_Ai(object):
""" L. Ai, K. Vafai, International Journal of Heat and Mass Transfer 49 (2006)
Table 2 Physiological parameters used in the numerical simulation
values given in milimeters
With D of endothelium depending on WSS"""

names       = [  'endothel' , 'intima', 'IEL'   ,'media' ]
D           = [   6e-11     , 5.0e-6  , 3.18e-9 , 5e-8   ]
V           = [   2.3e-2 ]*4
sigma       = [   0.9886    , 0.8292  , 0.8295  , 0.8660 ]
L           = [   2.        , 10.     ,    2.   , 200.   ]
k_react     = [   0.        , 0.      ,    0.   , 1.4e-4 ]
K           = [   3.2172e-15  , 2.2e-10   ,3.18e-13, 2e-12   ]
mu          = [   0.72e-3   , 0.72e-3 , 0.72e-3 , 0.72e-3 ]

def calculate_filration(self,dPmmHg):
mmHg2Pa = 133.3
dP = mmHg2Pa*dPmmHg
Rw = [L_*mu_/K_ for L_,K_,mu_ in zip(self.L,self.K,self.mu)]

self.Vfiltr = dP/sum(Rw)

def __init__(self,WSS=1.79):
""" Change units to mikrometers
Class can be initialized with a value of WSS in Pa
"""
dPmmHg=70

self.phi = phi(WSS=WSS)
self.sigma_end = sigma_end(phi=self.phi,w=14.3*1e-6,r_m = 11e-6)
self.Kend = Kend(w=14.3e-6,phi=self.phi)
self.K[0] = self.Kend*1e6
self.calculate_filration(dPmmHg=dPmmHg)

self.D = [D_*1e6 for D_ in self.D]
self.V = [ self.Vfiltr*1e6]*4
self.sigma[0] = self.sigma_end
self.D[0]=Diffusivity(w=14.3e-6,phi=self.phi,r_m = 11e-6)
def get_params(self):
return (self.D,self.V,self.sigma,self.L,self.k_react)

class LDL_Parameters_Vafai2006_Ai_without_D(object):
""" L. Ai, K. Vafai, International Journal of Heat and Mass Transfer 49 (2006)
Table 2 Physiological parameters used in the numerical simulation
values given in milimeters
With D of endothelium from that work.
"""
names       = [  'endothel' , 'intima', 'IEL'   ,'media' ]
D           = [   8.15e-11     , 5.0e-6  , 3.18e-9 , 5e-8   ]
V           = [   2.3e-2 ]*4
sigma       = [   0.9886    , 0.8292  , 0.8295  , 0.8660 ]
L           = [   2.        , 10.     ,    2.   , 200.   ]
k_react     = [   0.        , 0.      ,    0.   , 1.4e-4 ]
K           = [   3.2172e-15  , 2.2e-10   ,3.18e-13, 2e-12   ]
mu          = [   0.72e-3   , 0.72e-3 , 0.72e-3 , 0.72e-3 ]

def calculate_filration(self,dPmmHg):
mmHg2Pa = 133.3
dP = mmHg2Pa*dPmmHg
Rw = [L_*mu_/K_ for L_,K_,mu_ in zip(self.L,self.K,self.mu)]

self.Vfiltr = dP/sum(Rw)

def __init__(self,WSS=1.79):
""" Change units to mikrometers
Class can be initialized with a value of WSS in Pa
"""
dPmmHg=70

self.phi = phi(WSS=WSS)
self.sigma_end = sigma_end(phi=self.phi,w=14.3*1e-6,r_m = 11e-6)
self.Kend = Kend(w=14.3e-6,phi=self.phi)
self.K[0] = self.Kend*1e6
self.calculate_filration(dPmmHg=dPmmHg)

self.D = [D_*1e6 for D_ in self.D]
self.V = [ self.Vfiltr*1e6]*4
self.sigma[0] = self.sigma_end
def get_params(self):
return (self.D,self.V,self.sigma,self.L,self.k_react)

class LDL_Parameters_Olgac_WSS(object):
""" U. Olgac, V. Kurtcuoglu, D. Poulikakos, American Journal of Physiology-Heart and Circulatory
Physiology 294
The dependecy of WSS is implemented
"""
name       = [  'endothel' , 'wall' ]
D           = [   6e-11     ,   8.0e-7 ]
V           = [   2.3e-5 ]*2
sigma       = [   0.988    ,  0.8514   ]
L           = [   2.        , 338.  ]
k_react     = [   0.        ,  3.0e-4  ]
K           = [   3.32e-15  ,  1.2e-12]
mu          = [   0.72e-3   ,0.001]

def calculate_filration(self,dPmmHg):
mmHg2Pa = 133.3
dP = mmHg2Pa*dPmmHg
Rw = [L_*mu_/K_ for L_,K_,mu_ in zip(self.L,self.K,self.mu)]
self.Vfiltr = dP/sum(Rw)

def __init__(self,WSS=1.79):
""" Change units to mikrometers
Class can be initialized with a value of WSS in Pa
"""
dPmmHg=70

self.phi = phi(WSS=WSS)
self.sigma_end = sigma_end(phi=self.phi,w=14.3*1e-6,r_m = 11e-6)
self.Kend = Kend(w=14.3e-6,phi=self.phi)
self.K[0] = self.Kend*1e6
self.calculate_filration(dPmmHg=dPmmHg)

self.D = [D_*1e6 for D_ in self.D]
self.V = [ self.Vfiltr*1e6]*4
self.sigma[0] = self.sigma_end
self.D[0]=Diffusivity(w=14.3e-6,phi=self.phi,r_m = 11e-6)

def get_params(self):
return (self.D,self.V,self.sigma,self.L,self.k_react)



## Assembling and solving the discrete system¶

The LDL_Sim class is responsible for solve the differential equation.

In the first step the simulation region is discretized by the method discretize. In that function flux continuity between layers and boundary conditions are implemented. It is done in a following way:

• the layers are encoded as list layers in discretize, containing indices of boundaries for given discretization

• system parameters depending on space ($k,\vec u,\sigma,D_{eff}$) are sampled

• the diagonals of matrix are assembled using finite differences in space, neglecting for a moment the space dependence

• at regions boudaries the equation is replaced with $$J_L = J_R$$ it takes a form for a boundary at index j: $$c_j v_{j-1} - D_{j-1}\frac{c_j-c_{j-1}}{dx} = c_j v_{j+1} - D_{j+1}\frac{c_{j+1}-c_{j}}{dx}$$ collecting terms with the same $c$: $$c_j \left( v_{j-1}- v_{j+1} - \frac{D_{j-1}+D_{j+1}}{dx}\right) + c_{j-1}\frac{D_{j-1}}{dx} + c_{j+1}\frac{D_{j+1}}{dx} = 0$$

therefore we need to enforce the above equation on the system for each boundary inside the domain

• the system matrix is stored as scipy.sparse dia_matrix

The plot_c function is used for plot the concentration profiles.

In [4]:
class LDL_Sim(object):

def __init__(self, pars):
self.pars = pars
self.c_st = None
def discretize(self,N=2000):
self.N = N
k = np.ones(N)
v = np.ones(N)
Dyf = np.ones(N)

D,V,sigma,L,k_react = self.pars.get_params()
l = np.sum(L)
self.l = l
self.x=np.linspace(0,l,N)

layers=[0]+list( np.ceil( (N*(np.cumsum(L)/sum(L)))).astype(np.int32) )

for i,(l1,l2) in enumerate(zip(layers[:],layers[1:])):
k[l1:l2] = k_react[i]
v[l1:l2] = (1.0-sigma[i])*V[i]
Dyf[l1:l2] = D[i]
dx2_1 = (N-1)**2/l**2
dx_1 = (N-1)/l

diag_l = np.ones(N)*(np.roll(Dyf,-1)*dx2_1)
diag   = np.ones(N)*(-2.*Dyf*dx2_1 - k + v*dx_1)
diag_u = np.ones(N)*(np.roll(Dyf,1)*dx2_1 - np.roll(v,1)*dx_1)

# Layer's junctions
for j in layers[1:-1]:
diag[j] = v[j-1]-v[j+1]-(Dyf[j-1]+Dyf[j+1])*dx_1
diag_l[j-1] = Dyf[j-1]*dx_1
diag_u[j+1] = Dyf[j+1]*dx_1
#Boundary Conditions
diag[0] = 1
diag[-1] = 1
diag_u[0+1] = 0
diag_l[0-2] = 0

self.L = dia_matrix((np.array([diag_l,diag,diag_u]),np.array([-1,0,1])), shape=(N,N))

def solve_stationary(self,bc=[1,0]):
b = np.zeros(self.N)
b[0],b[-1] = bc
L = self.L.tocsr()
self.c_st = sp.sparse.linalg.linsolve.spsolve(L,b)

def plot_c(self,yrange=(0,0.2),xrange=(0,214),filename=None, color='red', alpha=0.2, style='-'):

i1,i2 = int(xrange[0]/self.l*self.N),int(xrange[1]/self.l*self.N)
plt.plot(self.x[i1:i2],self.c_st[i1:i2],color=color,linewidth=2, ls=style)
plt.ylim( *yrange)
plt.xlim( *xrange)
L=self.pars.L
d=[0]+np.cumsum(self.pars.L).tolist()
colors=['m','g','b','w']
for i,(l1,l2) in enumerate(zip(d[:],d[1:])):

plt.bar([l1,],yrange[1],l2-l1, color=colors[i], linewidth=0.3,  alpha=alpha)

plt.grid(True,axis='y', which='major')
plt.xlabel(r"$x \left[\mu m\right]$")
plt.ylabel(r"$c(x)$")
if filename!=None:
plt.savefig(filename)


## Class that performs simulation ¶

In that class the whole simulation is performed: parameters are initiated and then the equation is solved. As a result we get the object that contains the values of LDL concentration in discretized points and the function plot_c which can plot the concentration profile.

Class can be iitiated with name of one of the parameters sets:

• "4L_2012" for S. Chung, K. Vafai, International Journal of Biomechanics 45(2012)
• "4L_2006_Ai" for L. Ai, K. Vafai, International Journal of Heat and Mass Transfer 49 (2006) with D(WSS)
• "4L_2006_Ai_without_D" for L. Ai, K. Vafai, International Journal of Heat and Mass Transfer 49 (2006) with D from this publication
• "2L" for U. Olgac, V. Kurtcuoglu, D. Poulikakos, American Journal of Physiology-Heart and Circulatory Physiology 294
In [5]:
def LDL_simulation(wss=1.79, parameters="2012", bc=[1,0.0047],verbose=True):
if (parameters=="4L_2012"):
pars = LDL_Parameters_Vafai2012(WSS=wss)
elif (parameters=="4L_2006_Ai"):
pars = LDL_Parameters_Vafai2006_Ai(WSS=wss)
elif (parameters=="4L_2006_Ai_without_D"):
pars = LDL_Parameters_Vafai2006_Ai_without_D(WSS=wss)
elif (parameters=="2L"):
pars = LDL_Parameters_Olgac_WSS(WSS=wss)
else:
print "Parameters error"
return
sim = LDL_Sim(pars)
sim.discretize(130*214)
sim.solve_stationary(bc=bc)
if verbose:
print "The total surfaced LDL concentration:",np.sum(sim.c_st)*(sim.l/(sim.N-1))
return sim


# Examples¶

In this section the results presented in publications are calculated

## Small WSS=0.02¶

In [6]:
LDL_simulation(wss=0.02, parameters="4L_2012").plot_c(yrange=(0,5.0),xrange=(0,214), color='green', alpha=0.1)
LDL_simulation(wss=0.02, parameters="4L_2006_Ai").plot_c(yrange=(0,5.0),xrange=(0,214), color='blue', alpha=0.1, style='--')
legend(('4LC parametrs','4LA parametrs'));

The total surfaced LDL concentration: 74.2308568103
The total surfaced LDL concentration: 65.8469132801

In [7]:
LDL_simulation(wss=0.02, parameters="4L_2012").plot_c(yrange=(0,5.0),xrange=(0,25), color='green', alpha=0.1)
LDL_simulation(wss=0.02, parameters="4L_2006_Ai").plot_c(yrange=(0,5.0),xrange=(0,25), color='blue', alpha=0.1, style='--')
legend(('4LC parametrs','4LA parametrs'));

The total surfaced LDL concentration: 74.2308568103
The total surfaced LDL concentration: 65.8469132801


## High WSS=2.2¶

In [19]:
sim = LDL_simulation(wss=2.2, parameters="4L_2012")
sim.plot_c(yrange=(0,1.1),xrange=(0,214), color='green', alpha=0.1)
sim = LDL_simulation(wss=2.2, parameters="4L_2006_Ai")
sim.plot_c(yrange=(0,1.1),xrange=(0,214), color='blue', alpha=0.1, style='--')
legend(('4LC parametrs','4LA parametrs'));
# use e.g. xlim(0,10) to zoom

The total surfaced LDL concentration: 4.37625312868
The total surfaced LDL concentration: 4.28078869686

Out[19]:
<matplotlib.legend.Legend at 0x7fba91809f50>
In [16]:
sim = LDL_simulation(wss=2.2, parameters="4L_2012")
sim.plot_c(yrange=(0,1.1),xrange=(0,25), color='green', alpha=0.1)
sim = LDL_simulation(wss=2.2, parameters="4L_2006_Ai_without_D")
sim.plot_c(yrange=(0,1.1),xrange=(0,25), color='orange', alpha=0.1)
sim = LDL_simulation(wss=2.2, parameters="4L_2006_Ai")
sim.plot_c(yrange=(0,1.1),xrange=(0,25), color='blue', alpha=0.1, style='--')

legend(('4LC parametrs','4LA parametrs without D(WSS)','4LA parametrs with D(WSS)'));

The total surfaced LDL concentration: 4.37625312868
The total surfaced LDL concentration: 3.95620489143
The total surfaced LDL concentration: 4.28078869686


## Concentration in intima in function of WSS¶

In [10]:
WSSs = np.arange(0.0,3.0,0.1)
print WSSs
c_endo_wss=[LDL_simulation(wss=x, parameters="4L_2012",verbose=False).c_st[2*130] for x in WSSs]
c_endo_wss2=[LDL_simulation(wss=x, parameters="4L_2006_Ai",verbose=False).c_st[2*130] for x in WSSs]

[ 0.   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.   1.1  1.2  1.3  1.4
1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9]

In [11]:
plot (WSSs,c_endo_wss2)
plot (WSSs,c_endo_wss, ls='--')
title("WSS dependence of intima side LDL concentration at the endothelium", fontsize=10)
xlim([0.0,2.5])
xlabel('WSS [Pa]')
ylabel('$c_{w}end^{*}$')
grid(True)
legend(('4LA parametrs','4LC parametrs'))


Out[11]:
<matplotlib.legend.Legend at 0x7fba91979dd0>

# Two layer model ¶

## High WSS=2.2¶

In [12]:
sim = LDL_simulation(wss=2.2, parameters="2L")
sim.plot_c(yrange=(0,0.012),xrange=(0,340));

The total surfaced LDL concentration: 2.7040541751


## Low WSS=0.02¶

In [13]:
sim = LDL_simulation(wss=0.02, parameters="2L")
sim.plot_c(yrange=(0,0.25),xrange=(0,340))

The total surfaced LDL concentration: 14.9909714672


## Concentration in intima in function of WSS¶

In [14]:
c_endo_wss=[LDL_simulation(wss=x, parameters="2L",verbose=False).c_st[2*130] for x in WSSs]

In [20]:
plot (WSSs,c_endo_wss)
xlim([0.0,1.6])
ylim([0,0.3])
xlabel('WSS [Pa]')
ylabel('$c_{w}end^{*}$')
grid(True)
plt.title("WSS dependence of intima side LDL concentration at the endothelium", fontsize=10)

Out[20]:
<matplotlib.text.Text at 0x7fba918be5d0>