%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) 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 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) 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) 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 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')); 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)')); 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')) sim = LDL_simulation(wss=2.2, parameters="2L") sim.plot_c(yrange=(0,0.012),xrange=(0,340)); sim = LDL_simulation(wss=0.02, parameters="2L") sim.plot_c(yrange=(0,0.25),xrange=(0,340)) 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)