<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  WSS dependent parameters
The impact of WSS on the transport properties is given by the equations (8-11) developed by Olgac et al.
Step10:  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.
Step11:  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
Step12:  Class that performs simulation
In that class the whole simulation is performed
Step13:  Examples
In this section the results presented in publications are calculated
Small WSS=0.02
Step14:  High WSS=2.2
Step15:  Concentration in intima in function of WSS
Step16:  Two layer model
High WSS=2.2
Step17:  Low WSS=0.02
Step18:  Concentration in intima in function of WSS
<ASSISTANT_TASK:>
Python Code:
%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)
Explanation: 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.
End of explanation
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
Explanation: WSS dependent parameters
The impact of WSS on the transport properties is given by the equations (8-11) developed by Olgac et al.
End of explanation
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)
  
Explanation: 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.
End of explanation
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)
Explanation: 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.
End of explanation
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
Explanation: 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
End of explanation
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'));
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'));
Explanation: Examples
In this section the results presented in publications are calculated
Small WSS=0.02
End of explanation
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
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)'));
Explanation: High WSS=2.2
End of explanation
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]
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'))
 
Explanation: Concentration in intima in function of WSS
End of explanation
sim = LDL_simulation(wss=2.2, parameters="2L")
sim.plot_c(yrange=(0,0.012),xrange=(0,340));
Explanation: Two layer model
High WSS=2.2
End of explanation
sim = LDL_simulation(wss=0.02, parameters="2L")
sim.plot_c(yrange=(0,0.25),xrange=(0,340))
Explanation: Low WSS=0.02
End of explanation
c_endo_wss=[LDL_simulation(wss=x, parameters="2L",verbose=False).c_st[2*130] for x in WSSs]
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)
Explanation: Concentration in intima in function of WSS
End of explanation
<END_TASK>