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>