% The following section is to calculate the intracellular % concentration and permeability of each drug given pKa, logPn(o/c) and % electrical charges. % Clear the memory clear % Constant T = 373.15+37; % Body temperature (37centigrade) R = 8.314; % Universal gas constant F = 96484.56; % Faraday constant La = 0; % Lipid fraction in apical compartment Lc = 0.05; % Lipid fraction in cytosol Lm = 0; % Lipid fraction in mitochondria Lb = 0; % Lipid fraction in basolateral compartment Wa = 1-La; % Water fraction in apical compartment Wc = 1-Lc; % Water fraction in cytosol Wm = 1-Lm; % Water fraction in mitochondria Wb = 1-Lb; % Water fraction in basolateral compartment gamma_na = 1; % Activity coefficient of neutral molecules in apical compartment gamma_da = 1; % Activity coefficient of ionic molecules in apical compartment gamma_nc = 1.23026877; % Activity coefficient of neutral molecules in cytosol gamma_dc = 0.73799822; % Activity coefficient of ionic molecules in cytosol gamma_nm = 1; % Activity coefficient of neutral molecules in mitochondria gamma_dm = 1; % Activity coefficient of ionic molecules in mitochondria gamma_nb = 1; % Activity coefficient of neutral molecules in basolateral compartment gamma_db = 1; % Activity coefficient of ionic molecules in basolateral compartment Ca = 1 ; % Apical initical drug concentration (mM) % Areas and volumes (units in m^2 and m^3) Aa = 50*10^(-10) ; % The apical membrane surface area Aaa = 20*10^(-10) ; % The monolayer area Am = 100*3.14*10^(-12); % The mitochondrial membrane surface area Ab = 10^(-10); % The basolateral membrane surface area Vc = 10*10^(-15); % The cytosolic volume Vm = 100*5.24*10^(-19); % The mitochondrial volume Vb = 4.7*10^(-3); % The volume of basolateral compartment % Membrane potential (units in 'Voltage') Ea = -0.0093 ; % The membrane potential of apical membrane Em = -0.16; % The membrane potential of mitochondrial membrane Eb = 0.0119 ; % The membrane potential of basolateral membrane % pH values pHa = 6.8; % pH in apical compartment pHc = 7.0; % pH in cytosol pHm = 8.0; % pH in mitochondria pHb = 7.4; % pH in basolateral compartment % Read the drug properties % If the drug is neutral at physiological pH % the z is given 10^(-6) in stead of 0 since if z=0 the differenctial % equations can't be solved [DrugName,pKaall,logPnall,ZNall] = textread('drug.dat', '%s %f %f %f','commentstyle','matlab'); % Calculate the ionized logP(o/w); logPdall = logPnall-3.7 ; % The calculated results are saved in this file 'Peff_all.dat' len = length(pKaall) ; fid1 = fopen('Peff_all.dat','w'); str1 = ' Name --------------- pKa ----- logP_n,lip ---logP_d,lip---Cc(mM)-----Cm(mM)-------Cb(mM)-------Peff(cm/sec) ' ; fprintf(fid1,'%s\n',str1) ; for n = 1:len if ( abs(ZNall(n)-1) <= 10^(-6) ) logP_nlipT(n) = 0.33*logPnall(n)+2.2 ; logP_dlipT(n) = 0.37*logPdall(n)+2 ; end if ( abs(ZNall(n)+1) <= 10^(-6) ) logP_nlipT(n) = 0.37*logPnall(n)+2.2 ; logP_dlipT(n) = 0.33*logPdall(n)+2.6 ; end if ( abs(ZNall(n)-0) <= 10^(-5) ) logP_nlipT(n) = 0.33*logPnall(n)+2.2 ; logP_dlipT(n) = 0.33*logPdall(n)+2.2 ; end end % Get the first two decimals logP_nlip = round(logP_nlipT*100)/100 ; logP_dlip = round(logP_dlipT*100)/100 ; % Solve the differential equation system for each drug: % Given a system of linear ODE's expressed in matrix form: % Y' = AY+G with initial conditions Y(0) = RR, for n = 1:len pKa = pKaall(n); logP_n = logP_nlip(n) ; logP_d = logP_dlip(n) ; z = ZNall(n) ; % Parameters Calculation i = -sign(z) ; Na = ((z)*(Ea)*F)/(R*T); Nm = ((z)*(Em)*F)/(R*T); Nb = ((z)*(-Eb)*F)/(R*T); Pn = 10^(logP_n-6.7); Pd = 10^(logP_d-6.7); Kn_a = La*1.22*10^(logP_n); Kd_a = La*1.22*10^(logP_d); Kn_c = Lc*1.22*10^(logP_n); Kd_c = Lc*1.22*10^(logP_d); Kn_m = Lm*1.22*10^(logP_n); Kd_m = Lm*1.22*10^(logP_d); Kn_b = Lb*1.22*10^(logP_n); Kd_b = Lb*1.22*10^(logP_d); % Construct the matrix A and G fn_a = 1/(Wa/gamma_na+Kn_a/gamma_na+Wa*10^(i*(pHa-pKa))/gamma_da... +Kd_a*10^(i*(pHa-pKa))/gamma_da); fd_a = fn_a*10^(i*(pHa-pKa)); fn_c = 1/(Wc/gamma_nc+Kn_c/gamma_nc+Wc*10^(i*(pHc-pKa))/gamma_dc... +Kd_c*10^(i*(pHc-pKa))/gamma_dc); fd_c = fn_c*10^(i*(pHc-pKa)); fn_m = 1/(Wm/gamma_nm+Kn_m/gamma_nm+Wm*10^(i*(pHm-pKa))/gamma_dm... +Kd_m*10^(i*(pHm-pKa))/gamma_dm); fd_m = fn_m*10^(i*(pHm-pKa)); fn_b = 1/(Wb/gamma_nb+Kn_b/gamma_nb+Wb*10^(i*(pHb-pKa))/gamma_db... +Kd_b*10^(i*(pHb-pKa))/gamma_db); fd_b = fn_b*10^(i*(pHb-pKa)); k11 = -(Aa/Vc)*Pn*fn_c-(Aa/Vc)*Pd*Na*fd_c*exp(Na)/(exp(Na)-1)... -(Am/Vc)*Pn*fn_c-(Am/Vc)*Pd*Nm*fd_c/(exp(Nm)-1)... -(Ab/Vc)*Pn*fn_c-(Ab/Vc)*Pd*Nb*fd_c/(exp(Nb)-1) ; k12 = (Am/Vc)*Pn*fn_m+(Am/Vc)*Pd*Nm*fd_m*exp(Nm)/(exp(Nm)-1) ; k13 = (Ab/Vc)*Pn*fn_b+(Ab/Vc)*Pd*Nb*fd_b*exp(Nb)/(exp(Nb)-1) ; S1 = (Aa/Vc)*Ca*(Pn*fn_a+Pd*Na*fd_a/(exp(Na)-1)) ; k21 = (Am/Vm)*Pn*fn_c+(Am/Vm)*Pd*Nm*fd_c/(exp(Nm)-1) ; k22 = -(Am/Vm)*Pn*fn_m-(Am/Vm)*Pd*Nm*fd_m*exp(Nm)/(exp(Nm)-1) ; k23 = 0; S2 = 0; k31 = (Ab/Vb)*Pn*fn_c+(Ab/Vb)*Pd*Nb*fd_c/(exp(Nb)-1) ; k32 = 0; k33 = -(Ab/Vb)*Pn*fn_b-(Ab/Vb)*Pd*Nb*fd_b*exp(Nb)/(exp(Nb)-1) ; S3 = 0; A = [k11, k12, k13; k21, k22, k23; k31, k32, k33]; G = [S1, S2, S3]'; RR = [0,0,0]'; t = 1000; % Calculate the intracellular concentration and permeability and t=1000s % which is at steady state [V,E] = eig(A); E = diag(E); H = inv(V)*G; B = V \ RR; C = B + H./E; Z = -(H./E) + exp(t * E).*C ; Y = real(V * Z); Y = Y'; Peff = Y(3)*Vb/(t*Aaa*Ca); NA = [pKa, logP_n, logP_d, Y, Peff*10^(8)]; str = DrugName{n}; fprintf(fid1,'%s\t %12.2f %12.2f %12.2f %12.2f %12.2f %+12.4e %12.2f\n',str, NA') ; end fclose(fid1);