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)
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
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)
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 [ ]:
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)
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
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'));
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'));
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
Out[19]:
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)'));
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]
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]:
In [12]:
sim = LDL_simulation(wss=2.2, parameters="2L")
sim.plot_c(yrange=(0,0.012),xrange=(0,340));
In [13]:
sim = LDL_simulation(wss=0.02, parameters="2L")
sim.plot_c(yrange=(0,0.25),xrange=(0,340))
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]:
In [15]: