% The following section is to calculate the intracellular concentration and % permeability of an arbitary chemical space % Change logPn and logPd independently % 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 % Parameters Calculation i = -sign(z) ; Na = ((z)*(Ea)*F)/(R*T); Nm = ((z)*(Em)*F)/(R*T); Nb = ((z)*(-Eb)*F)/(R*T); N = zeros(1,7); for pKa = 1:0.15:14 for logP_n = -5:0.15:5 for logP_d = -5:0.15:5 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); 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; % Solve the system. [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); N = [N; pKa, logP_n, logP_d, Y, Peff]; end end end N = N(2:length(N),:); % output results into file 'Total_Chem.dat' fid = fopen('Total_Chem.dat','w'); str = ' pKa ------------ logP_n ------------ logP_d ------------ Cc ------------ Cm ------------ Cb ------------ Peff' ; fprintf(fid,'%s\n',str) ; fprintf(fid,'%+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e\n',N') ; fclose(fid); % The following section is to construct the Permeant-Nontoxic chemical % space based on the above total chemical space % Clear the memory clear % Calculate the intracellular concentration and permeability of metoprolol % (the reference drug). % 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 % Metoprolol properties z = 1 ; pKa = 9.7 ; logPn = 1.88 ; logPd = logPn-3.7 ; logP_n = round((0.33*logPn+2.2)*100)/100 ; logP_d = round((0.37*logPd+2)*100)/100 ; % 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); 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; [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); % Storage of calculated intracellular concentration and permeability of metoprolol NA = [Y, Peff]; % Construct the Permeant-Nontoxic chemical space by comparing with % metoprolol % Read the data filename = 'Total_Chem.dat' ; fid = fopen(filename) ; fgetl(fid) ; dataall = fscanf(fid,'%f %f %f %f %f %f %f',[7,inf]) ; dataall = dataall' ; datapts = dataall(1:length(dataall),1:7) ; L = length(datapts); NN1 = zeros(1,7); NN2 = zeros(1,7); % definition of permeant-nontoxic chemical space % permeability is equal or larger than metoprolol permeability % cytosolic and mitochondrial concentration are equal or less than % metoprolol cytosolic and mitochondrial concentration for i = 1:L if ( datapts(i,7) >= NA(1,4) && datapts(i,4) <= NA(1,1) && datapts(i,5) <= NA(1,2) ) NN1 = [NN1; datapts(i,:)]; else NN2 = [NN2; datapts(i,:)]; end end NN1 = NN1(2:length(NN1),:); NN2 = NN2(2:length(NN2),:); % output permeant-nontoxic results into file 'PermNontox.dat' fid1 = fopen('PermNontox.dat','w'); str1 = ' pKa ------------ logP_n ------------ logP_d ------------ Cc ------------ Cm ------------ Cb ------------ Peff' ; fprintf(fid1,'%s\n',str1) ; fprintf(fid1,'%+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e\n',NN1') ; fclose(fid1); % output permeant-toxic results into file 'PermTox.dat' fid1 = fopen('PermTox.dat','w'); str1 = ' pKa ------------ logP_n ------------ logP_d ------------ Cc ------------ Cm ------------ Cb ------------ Peff' ; fprintf(fid1,'%s\n',str1) ; fprintf(fid1,'%+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e %+12.8e\n',NN2') ; fclose(fid1); % The following section is to minimize data for plotting % clear memories ; clear ; % read data ; filename = 'PermNontox.dat' ; % file's name fid = fopen(filename) ; % open file fgetl(fid) ; % read an useless line: the first line of the file is useless dataall = fscanf(fid,'%f %f %f %f %f %f %f',[7,inf]) ; dataall = dataall' ; % dataall' is the transpose of dataall datapts = dataall(1:length(dataall),1:3) ; % only the first 3 columns will be used datapts = sortrows(datapts) ; [NR, NC] = size(datapts) ; INCOM = zeros(1,3); %Storage for incomplete data for n = 1:0.15:14 % pka TEMP1 = zeros(1,3) ; for i = 1:NR if (abs(datapts(i,1)-n) < 1e-6 ) TEMP1 = [TEMP1; datapts(i,:)]; end end [TNR1, TNC1] = size(TEMP1) ; if (TNR1 > 1 ) TEMP1 = TEMP1(2:TNR1,:) ; for k = -5:0.15:5 %logPn TEMP2 = zeros(1,3) ; for l = 1:size(TEMP1) if (abs(TEMP1(l,2)-k) < 1e-6 ) TEMP2 = [TEMP2; TEMP1(l,:)]; end end [TNR, TNC] = size(TEMP2) ; if (TNR >= 3) INCOM = [INCOM; TEMP2(2,:); TEMP2(TNR,:)] ; end if ( abs(TNR-2) < 1e-6 ) INCOM = [INCOM; TEMP2(2,:)] ; end if ( abs(TNR-1) < 1e-6 ) INCOM = [INCOM] ; end end end end INCOM = INCOM(2:size(INCOM),:); % output results into file fid1 = fopen('PermNontox_incomp.dat','w'); str1 = ' pKa ------------ logP_n ------------ logP_d ' ; fprintf(fid1,'%s\n',str1) ; fprintf(fid1,'%+12.8e %+12.8e %+12.8e\n',INCOM') ; fclose(fid1); % The following section is to plot permeant-nontoxic chemical spaces % clear memories ; clear ; % read data ; filename = 'PermNontox_incomp.dat' ; % file's name fid = fopen(filename) ; % open file fgetl(fid) ; % read an useless line: the first line of the file is useless dataall = fscanf(fid,'%f %f %f',[3,inf]) ; dataall = dataall' ; % dataall' is the transpose of dataall datapts = dataall(1:length(dataall),1:3) ; % only the first 3 columns will be used % analyse the data num_x = 0 ; % number of planes with the same x (pKa) num_y = zeros(1000,1) ; % number of lines with the same y (logP_n) [on one plane with the same x (pKa)] xyz_up = zeros(1000,1000,3) ; % the coordinates of the upper points. xyz_up(m,n,k) corresponds to the upper % point on the m-th plane [with the same x % (pKa)] and the n-th line [with the same y % (logP_n)]. k=1 means the x coordinate, % k=2 means the y coordinate, k=3 means the % z coordinate (logP_d). xyz_dn = zeros(1000,1000,3) ; % the coordinates of the lower points. index_up = zeros(1000,1000) ; % index of the upper points index_dn = zeros(1000,1000) ; % index of the lower points figure(1) ; clf ; % initialize the figure hold on ; xlabel('pK_a') ; ylabel('logP_n') ; zlabel('logP_d') ; set(gca,'XLim',[0 15],'YLim',[-5 5],'ZLim',[-5 5]); title(' Permeant Nontoxic Chemical Space ') ; set(gca,'LineWidth',1.0,'FontSize',20,'FontWeight','Demi') ; set(get(gca,'XLabel'),'FontSize',20,'FontWeight','Demi') ; set(get(gca,'YLabel'),'FontSize',20,'FontWeight','Demi') ; set(get(gca,'ZLabel'),'FontSize',20,'FontWeight','Demi') ; set(get(gca,'Title'),'FontSize',20,'FontWeight','Demi') ; flag_x = 0 ; flag_y = 0 ; z_up_value = -inf ; z_dn_value = inf ; len = length(datapts) ; ind = 0 ; x_value = inf ; while( ind < len) ind = ind + 1; pt = datapts(ind,:) ; % the ind-th point if (x_value ~= pt(1)) % select the new plane with the same x (pKa) num_x = num_x+1 ; x_value = pt(1) ; y_value = inf ; end if (y_value ~= pt(2)) % select the new line with the same y (logP_n) [in the plane with the same x (pKa)] num_y(num_x) = num_y(num_x) + 1 ; y_value = pt(2) ; z_up_value = -inf ; z_dn_value = inf ; end % set the x and y coordinates for the upper and lower points xyz_up(num_x,num_y(num_x),1) = pt(1) ; xyz_dn(num_x,num_y(num_x),1) = pt(1) ; xyz_up(num_x,num_y(num_x),2) = pt(2) ; xyz_dn(num_x,num_y(num_x),2) = pt(2) ; % set the z coordinates for the upper and lower points if (pt(3) > z_up_value) z_up_value = pt(3) ; xyz_up(num_x,num_y(num_x),3) = pt(3) ; index_up(num_x,num_y(num_x)) = ind ; end if (pt(3) < z_dn_value) z_dn_value = pt(3) ; xyz_dn(num_x,num_y(num_x),3) = pt(3) ; index_dn(num_x,num_y(num_x)) = ind ; end end for m = 1:num_x % the m-th plane n = 1 ; x1 = [xyz_up(m,n,1) xyz_dn(m,n,1)] ; y1 = [xyz_up(m,n,2) xyz_dn(m,n,2)] ; z1 = [xyz_up(m,n,3) xyz_dn(m,n,3)] ; plot3(x1,y1,z1) ; % plot the first line which connects the upper and the lower points n = num_y(m) ; x1 = [xyz_up(m,n,1) xyz_dn(m,n,1)] ; y1 = [xyz_up(m,n,2) xyz_dn(m,n,2)] ; z1 = [xyz_up(m,n,3) xyz_dn(m,n,3)] ; plot3(x1,y1,z1) ; % plot the last line which connects the upper and the lower points end m = 1 ; for n = 1:num_y(m)-1 x1 = [xyz_up(m,n,1) xyz_up(m,n+1,1)] ; y1 = [xyz_up(m,n,2) xyz_up(m,n+1,2)] ; z1 = [xyz_up(m,n,3) xyz_up(m,n+1,3)] ; plot3(x1,y1,z1) ; % plot the line which connects the upper points only x1 = [xyz_dn(m,n,1) xyz_dn(m,n+1,1)] ; y1 = [xyz_dn(m,n,2) xyz_dn(m,n+1,2)] ; z1 = [xyz_dn(m,n,3) xyz_dn(m,n+1,3)] ; plot3(x1,y1,z1) ; % plot the line which connects the lower points only end m = num_x ; for n = 1:num_y(m)-1 x1 = [xyz_up(m,n,1) xyz_up(m,n+1,1)] ; y1 = [xyz_up(m,n,2) xyz_up(m,n+1,2)] ; z1 = [xyz_up(m,n,3) xyz_up(m,n+1,3)] ; plot3(x1,y1,z1) ; % plot the line which connects the upper points only x1 = [xyz_dn(m,n,1) xyz_dn(m,n+1,1)] ; y1 = [xyz_dn(m,n,2) xyz_dn(m,n+1,2)] ; z1 = [xyz_dn(m,n,3) xyz_dn(m,n+1,3)] ; plot3(x1,y1,z1) ; % plot the line which connects the lower points only end for m = 1:num_x-1 % plot the outer lines connects the upper (lower) points with different x x1 = [xyz_up(m,1,1) xyz_up(m+1,1,1)] ; y1 = [xyz_up(m,1,2) xyz_up(m+1,1,2)] ; z1 = [xyz_up(m,1,3) xyz_up(m+1,1,3)] ; plot3(x1,y1,z1) ; x1 = [xyz_up(m,num_y(m),1) xyz_up(m+1,num_y(m+1),1)] ; y1 = [xyz_up(m,num_y(m),2) xyz_up(m+1,num_y(m+1),2)] ; z1 = [xyz_up(m,num_y(m),3) xyz_up(m+1,num_y(m+1),3)] ; plot3(x1,y1,z1) ; x1 = [xyz_dn(m,1,1) xyz_dn(m+1,1,1)] ; y1 = [xyz_dn(m,1,2) xyz_dn(m+1,1,2)] ; z1 = [xyz_dn(m,1,3) xyz_dn(m+1,1,3)] ; plot3(x1,y1,z1) ; x1 = [xyz_dn(m,num_y(m),1) xyz_dn(m+1,num_y(m+1),1)] ; y1 = [xyz_dn(m,num_y(m),2) xyz_dn(m+1,num_y(m+1),2)] ; z1 = [xyz_dn(m,num_y(m),3) xyz_dn(m+1,num_y(m+1),3)] ; plot3(x1,y1,z1) ; end figure(2) ; clf ; % initialize the figure hold on ; xlabel('pK_a') ; ylabel('logP_n') ; zlabel('logP_d') ; set(gca,'XLim',[0 15],'YLim',[-5 5],'ZLim',[-5 5]); set(gca,'LineWidth',5.0,'FontSize',20,'FontWeight','Bold','FontName','Times') ; set(get(gca,'XLabel'),'FontSize',20,'FontWeight','Bold','FontName','Times') ; set(get(gca,'YLabel'),'FontSize',20,'FontWeight','Bold','FontName','Times') ; set(get(gca,'ZLabel'),'FontSize',20,'FontWeight','Bold','FontName','Times') ; FaceVertexCDataValue = 0.7 ; FaceAlphaValue = 0.6 ; % between 0-1, 0 means transparent, 1 opaque for m = 1:num_x-1 % plot the outer lines connects the upper (lower) points with different x % the outer one with the minimal y index = zeros(1,4) ; index(1) = index_dn(m,1) ; index(2) = index_dn(m+1,1) ; index(3) = index_up(m+1,1) ; index(4) = index_up(m,1) ; patch('Faces',index,'Vertices',datapts,'FaceVertexCData',FaceVertexCDataValue,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ; % the outer one with the maximal y index = zeros(1,4) ; index(1) = index_dn(m,num_y(m)) ; index(2) = index_dn(m+1,num_y(m+1)) ; index(3) = index_up(m+1,num_y(m+1)) ; index(4) = index_up(m,num_y(m)) ; patch('Faces',index,'Vertices',datapts,'FaceVertexCData',FaceVertexCDataValue,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ; % the lower one with the minimal z index = zeros(1,num_y(m)+num_y(m+1)) ; for n = 1:num_y(m) index(n) = index_dn(m,n) ; end for n = 1:num_y(m+1) index(num_y(m)+num_y(m+1)+1-n) = index_dn(m+1,n) ; end patch('Faces',index,'Vertices',datapts,'FaceVertexCData',FaceVertexCDataValue,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ; % the upper one with the maximal z index = zeros(1,num_y(m)+num_y(m+1)) ; for n = 1:num_y(m) index(n) = index_up(m,n) ; end for n = 1:num_y(m+1) index(num_y(m)+num_y(m+1)+1-n) = index_up(m+1,n) ; end patch('Faces',index,'Vertices',datapts,'FaceVertexCData',FaceVertexCDataValue,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ; end m = 1 ; % for the plane with the minimal x index = zeros(1,2*num_y(m)) ; for n = 1:num_y(m) index(n) = index_dn(m,n) ; index(2*num_y(m)+1-n) = index_up(m,n) ; end patch('Faces',index,'Vertices',datapts,'FaceVertexCData',FaceVertexCDataValue,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ; m = num_x ; % for the plane with the maximal x index = zeros(1,2*num_y(m)) ; for n = 1:num_y(m) index(n) = index_dn(m,n) ; index(2*num_y(m)+1-n) = index_up(m,n) ; end patch('Faces',index,'Vertices',datapts,'FaceVertexCData',0.7,... 'FaceColor',[190,190,190]/256,'FaceAlpha',FaceAlphaValue) ;