eeg_lap_test_script - explore the vector calculus of scalp current density. It starts with a potential field (ie, scalp voltage), which is a scalar field in x,y. We then calculate the first directional derivative, called the gradient. This is a vector field with unit magnitude and direction of greatest incline of the potential field (ie, all vectors are directed away from a negative potential peak and toward a positive potential peak). The next step is to calculate the current density vector. This depends on the conductivity (sigma) of the material from which the potential field is measured, given by: current density = -sigma * grad(v) Lets assume the scalp potential is measured from a scalp with constant conductivity and a good estimate of this value is 0.35 Siemens/m (0.035 Siemens/cm). Given a gradient of the potential field over meters, we have the following: scalp current density = -0.35 * grad(v) The next step is to calculate the divergence of the gradient vector. This is a scalar field that indicates the amount of expansion (+ve) or contraction (-ve) in the direction of the gradient field. div( scalp current density ) = -0.35 * div( grad(v) ) Mathematically, this step is performed with the Laplacian operator (del^2). The Laplacian is equivalent to the divergence of the gradient - expressed as: div( grad(v) ) or del . del(v) or del^2(v) However, for a Guassian scalp surface, with no current sources/sinks, the divergence of the gradient is zero, because no current can be conducted along a scalp normal direction from the scalp into the air. To solve this problem we need the sum of the partial derivates of scalp potential in x and y (tangential components). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 % eeg_lap_test_script - explore the vector calculus of scalp current density. 0002 % 0003 % It starts with a potential field (ie, scalp voltage), 0004 % which is a scalar field in x,y. 0005 % 0006 % We then calculate the first directional derivative, called 0007 % the gradient. This is a vector field with unit magnitude and 0008 % direction of greatest incline of the potential field (ie, 0009 % all vectors are directed away from a negative potential peak 0010 % and toward a positive potential peak). 0011 % 0012 % The next step is to calculate the current density vector. This 0013 % depends on the conductivity (sigma) of the material from which the 0014 % potential field is measured, given by: 0015 % 0016 % current density = -sigma * grad(v) 0017 % 0018 % Lets assume the scalp potential is measured from a scalp 0019 % with constant conductivity and a good estimate of this value 0020 % is 0.35 Siemens/m (0.035 Siemens/cm). Given a gradient of the 0021 % potential field over meters, we have the following: 0022 % 0023 % scalp current density = -0.35 * grad(v) 0024 % 0025 % The next step is to calculate the divergence of the gradient 0026 % vector. This is a scalar field that indicates the amount of 0027 % expansion (+ve) or contraction (-ve) in the direction of the 0028 % gradient field. 0029 % 0030 % div( scalp current density ) = -0.35 * div( grad(v) ) 0031 % 0032 % Mathematically, this step is performed with the Laplacian 0033 % operator (del^2). The Laplacian is equivalent to the divergence 0034 % of the gradient - expressed as: 0035 % 0036 % div( grad(v) ) or del . del(v) or del^2(v) 0037 % 0038 % However, for a Guassian scalp surface, with no current sources/sinks, 0039 % the divergence of the gradient is zero, because no current can 0040 % be conducted along a scalp normal direction from the scalp into the 0041 % air. 0042 % 0043 % To solve this problem we need the sum of the partial derivates of 0044 % scalp potential in x and y (tangential components). 0045 % 0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0047 0048 clear; close all; 0049 0050 SHADE = 'interp'; % FLAT, INTERP, {FACETED} 0051 VIEW = [-20,20]; 0052 0053 rot = 0; % Zaxis label rotation 0054 align = 'right'; 0055 0056 fig1 = figure('Position',[10 10 1200 900]), rotate3d; 0057 sub11 = subplot(3,2,1); sub12 = subplot(3,2,2); sub13 = subplot(3,2,3); 0058 sub14 = subplot(3,2,4); sub15 = subplot(3,2,5); sub16 = subplot(3,2,6); 0059 0060 fig2 = figure('Position',[500 100 500 800]); 0061 sub21 = subplot(2,1,1); sub22 = subplot(2,1,2); 0062 0063 % A scalar potential field in x,y is any function f(x,y), eg: 0064 0065 MIN = -pi/2; MAX = pi/2; grid = pi/10; 0066 [x,y] = meshgrid(MIN:grid:MAX); r = x.^2 + y.^2; 0067 0068 % Model a charge distributions 0069 Q1 = sin(x); Q2 = sin(y); z = (Q1 .* Q2); 0070 % cubic spline interpolation of charge density rho(x,y) 0071 zi = interp2(z,1,'cubic'); xi = linspace(MIN,MAX,size(zi,1)); yi = xi; 0072 0073 figure(fig1); subplot(sub11); surfc(xi,yi,zi), shading (SHADE), hold on, colorbar, view(VIEW); 0074 title 'Charge Density (C/m^2)', zlabel('coulombs (C)','Rotation',rot,'HorizontalAlignment',align), axis tight 0075 0076 % Electric Field Strength 0077 0078 %e0 = 8.854E-12; % permittivity of free space 0079 %k = 1 / (4 * pi * e0); % permittivity constant (8.9877e+009) 0080 0081 [u,v,en] = surfnorm(x,y,z); % 3D surface normals 0082 en = en.*z.*r; % orient normal in direction of charge & scale by distance 0083 0084 [ex,ey] = gradient(-z,grid,grid); % electric field, gradient of -charge (del[-z]) 0085 ex = ex.*r; % scale ex component by distance 0086 ey = ey.*r; % scale ey component by distance 0087 0088 figure(fig1); subplot(sub12); surfc(xi,yi,zi), shading (SHADE), hold on, colorbar, view(VIEW) 0089 title 'Electric Field (E, newtons/coulomb * 8.9877e09)', zlabel('N/C','Rotation',rot,'HorizontalAlignment',align) 0090 quiver3(x,y,z,ex,ey,en,0), axis tight 0091 0092 figure(fig2); subplot(sub21); contour(xi,yi,zi), hold on, quiver(x,y,ex,ey,0), hold off 0093 title 'Electric Field Vectors (N/C)' 0094 0095 % Electric potential (volts = joules/coulomb) is proportional to minus Electric Field Strength 0096 % and equivalent to the charge density in free space 0097 0098 figure(fig1); subplot(sub13); surfc(xi,yi,zi), shading (SHADE), hold on, colorbar, view(VIEW); 0099 title 'Electric Potential (V; 1 volt = 1 joule/coulomb)', axis tight 0100 zlabel('Volt','Rotation',rot,'HorizontalAlignment',align) 0101 0102 0103 % Current Density = Conductivity * Electric Field Strength 0104 0105 % For a material conducting medium, current density is equal to 0106 % electric field strength multiplied by conductivity (sigma = 1/resistivity). 0107 0108 sigma = 0.35; % conductivity of human scalp tissues 0109 0110 cdx = ex.*sigma; % scale ex component by conductivity 0111 cdy = ey.*sigma; % scale ey component by conductivity 0112 cdn = en.*sigma; % scale en component by conductivity 0113 0114 figure(fig1); subplot(sub14); surfc(xi,yi,zi), shading (SHADE), hold on, colorbar, view(VIEW) 0115 title 'Current Density (Jscalp = E * 0.35)', zlabel('Amp/m^2','Rotation',rot,'HorizontalAlignment',align) 0116 quiver3(x,y,z,cdx,cdy,cdn,0), axis tight 0117 0118 figure(fig2); subplot(sub22); contour(xi,yi,zi), hold on, quiver(x,y,cdx,cdy,0), hold off 0119 title 'Current Density 2D Gradient' 0120 0121 % The Laplacian of a scalar field, eg electric potential (voltage) 0122 0123 lap = 4*del2(z,grid,grid); % Hjorth, nearest neighbour Laplacian 0124 lapi = interp2(lap,1,'cubic'); 0125 0126 figure(fig1); subplot(sub15); surfc(xi,yi,lapi), shading (SHADE), colorbar, axis tight, view(VIEW) 0127 title 'Laplacian of Electric Potential (Lap(V) = div[grad[V]])' 0128 zlabel('Volt/m^2','Rotation',rot,'HorizontalAlignment',align); 0129 0130 %figure('Name','Potential, Coloured by Laplacian','Position',[650 10 600 400]) 0131 %surfc(xi,yi,zi,lapi), shading (SHADE), colorbar, axis tight, view(VIEW) 0132 0133 % The divergence of the current density; Div(Current Density) = - conductivity * Laplacian(voltage) 0134 0135 scdi = -0.35 * lapi; 0136 0137 figure(fig1); subplot(sub16); surfc(xi,yi,scdi), shading (SHADE), colorbar, axis tight, view(VIEW) 0138 title 'Divergence of Current Density (div(J) = -0.35 * Lap(V))', 0139 zlabel('Amp/m^2','Rotation',rot,'HorizontalAlignment',align);