电极优化方法梳理
电场建模
考虑图1中具有标量电导率场$σ$的非均匀体积（我们忽略各向异性，并假设表1中列出的各向同性电导率）。由于组织没有净电流源或汇，因此组织内的电流密度$J$具有零发散：$▽·J=0(▽=[\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}] )$。因此，当通过电极向该体积的边界施加电流时，体积中产生的电势分布$V$可作为拉普拉斯方程的解：
假定当在组织边界的电导率变化时，电位和电流密度是连续的。通常，该解不具有闭合形式。然而，可以通过将体积离散为一组有限元（每个有限元具有单个电导率值）并使用有限元法（FEM）求解（1）来找到数值近似。假设使用$M$个电极和一个额外的参考电极，将$B_m，m=1,…,M$定义为第$m$组边界条件。利用$FEM$可以得出对于每个边界条件$B_m$的有效电阻率$a_m∈R^3$，将正负电极之间的电流流动和贯穿体积的场强$e_m∈R^3$相联系：
我们使用术语“有效电阻率”来表示将施加电流的大小与感应电场相关的量。该有效电阻率当然与组成组织的各个电阻率（或电导率）有关
其中$e_m(r_n)$是在$FEM$节点$r_n∈R^3, n=1,…,N$处由电极$m$刺激产生的电场向量；$s_m$是施加的电流密度幅值（为简单起见，我们假设施加电流的方向是固定的，例如，垂直于电极表面，并以标量幅值$s_m$工作），$N$是$FEM$节点的个数。
我们现在同时刺激所有电极，使得电极$m$处的电流密度大小由$s_m$给出，然后拉普拉斯方程的线性特点决定位置$r_n$处的净电场由下式给出：
在离散位置空间中按列堆叠并以矩阵形式重新写成：
其中：
净电场是由$FEM$求解器计算的每个双极配置产生的单个电场的线性组合。显然，场的性质与单个电流密度$s_m$密切相关。这一事实使我们能够调整这些电流密度，以使产生的电场相对于指定的测量值是最佳的。选择系数$s_m$来塑造感应场的问题类似于阵列信号处理中的“波束形成”问题。
% 提取A矩阵，以及得到总位置点数Nlocs(roast_target.m line 498-515)
% extract A matrix corresponding to the brain
indBrain = elem((elem(:,5)==1 | elem(:,5)==2),1:4);
indBrain = unique(indBrain(:));
A = A_all(indBrain,:,:);
% convert pseudo-world coordinates back to voxel coordinates for targeting,
% as targeting code works in the voxel space
nodeV = zeros(size(node,1),3);
for i=1:3, nodeV(:,i) = node(:,i)/hdrInfo.pixdim(i); end
locs = nodeV(indBrain,1:3);
isNaNinA = isnan(sum(sum(A,3),2)); % make sure no NaN is in matrix A or in locs
if any(isNaNinA), A = A(~isNaNinA,:,:); locs = locs(~isNaNinA,:); end
Nlocs = size(locs,1);
p.Nlocs = Nlocs;
Nelec = size(A,3);
A = reshape(A,Nlocs*3,Nelec);
电极优化
Weight Least squares（wls）
参考文献：Optimized multi-electrode stimulation increases focality and intensity at target
所需场强定义为$e_d(r_n)∈R^3$：
其中$e_0$定义为目标区域的期望方向和强度。$Τ$是定义目标区域的一系列节点。$Τ^c$是$\Tau$的补集。$e_o$指定方向的合理选择是与颅骨表面径向或切向，特别是在皮质刺激的情况下。
最小化期望场与（4）的线性叠加可实现的场之间的平方误差，可得到经典的最小二乘解：
定义对角矩阵$W∈R^{3N×3N}$，其中非零元素设置为：
其中$[·]$将参数舍入为最接近的整数。最终的约束优化问题可以写为：
% 最小二乘法程序实现（in ls_l1.m）
function [x,cvx_status] = ls_l1(A,d,ub,verbose)
% A：A矩阵
% d：目标场强
% ub: 最大允许电流
% verbose：显示进度详细信息
n = size(A,2);
if verbose
cvx_begin
else
cvx_begin quiet
end
variable x(n);
minimize( norm(A*x-d,2) );
subject to
norm([x;-sum(x)],1) <= 2*ub; cvx_end < code>% w权重矩阵的获得（in optimize_prepare.m）
targetRadius = p.targetRadius; % 目标半径，默认为2
target_nodes = cell(numOfTargets,1);
for n = 1:numOfTargets
target_nodes{n} = find(distances_to_target(:,n)% wls，调用ls_l1.m (in optimize_currents.m)
% [U,S,V] = svd(repmat(sqrt(w),1,size(A,2)).*A, 0);（optimize_prepare.m line 80）
% svd奇异值分解：https://zhuanlan.zhihu.com/p/29846048
sqrtw = sqrt(w); % w为权重矩阵
[s_opt,status] = ls_l1(S*V',U'*(sqrtw.*d),S_max,verbose);
x_opt = A*s_opt;
Weighted least squares with individual $l_1$ constraint （wls-per）
限制总电流的安全约束由注入电流的效应是全局相加的概念下得出的。然而，这可能不是电刺激安全限制的准确表示，因为某些电极配置可能仅在局部起作用；例如，头部一侧的双极电极不会影响头部另一侧的脑组织。在这种情况下，由每个电极正下方的皮肤反应和感觉水平引导的更局部的安全约束可能更合适。在这些情况下，放松$||s||_1$上的安全约束，并且限制每个单独电极上的电流绝对值可能是合适的。如果电极相距相当大的距离，这种程序可能特别有效。相应的优化问题写为：
数值计算将显示，产生的电场比最小二乘法产生的电场更强烈，但代价是流过大脑的电流更大。
% 添加单个电极电流约束的加权最小二乘法wls-per.m (in optimize_currents.m)
function [x,cvx_status] = ls_l1per(A,d,ub,verbose)
n = size(A,2);
if verbose
cvx_begin
else
cvx_begin quiet
end
variable x(n);
minimize( norm(A*x-d,2) );
subject to
norm([x;-sum(x)],1) <= 2="1mA" 2*ub; norm([x;-sum(x)],inf) <="ub/2;" % 单个电极的电流约束，ub="2mA," ub cvx_end code>Linearly constrained minimum variance （LCMV）
在波束形成文献中，一种常见的方案是在利用剩余自由度最小化总功率的同时强制执行硬线性约束（例如，强制给定方向的期望增益）。类似地，我们希望在单个节点上获得指定的电场，同时最小化其他节点的电场。这样的过程不需要任何加权，并且硬约束确保我们在目标处获得所需的场。从满足硬约束的所有电流分布$s$中，线性约束最小方差（LCMV）程序选择整个体积中具有最低总电场功率的电流分布$||As||^2$：
$n_0$是目标节点参数。通过拉格朗日乘子法可以得出：
然而，请注意，不能保证在目标处产生的电场达到最大值；相比之下，LCMV的一个已知缺点是在非目标区域形成不期望的“旁瓣”。这是为确保目标处准确的所需场强和方向而付出的代价。此外，如果选择目标节点处的期望强度太大，则所需电流可能违反安全标准；在这种情况下，必须为目标强度选择较低的值。
为确保刺激的安全性，添加总$l_1$约束以及个体$l_1$约束的LCMV优化问题为：
% 添加总l1约束以及个体l1约束的LCMV算法实现（lcmv_per.m）
function [x,cvx_status] = lcmv_l1per(A,C,f,ub,verbose)
n = size(A,2);
if verbose
cvx_begin
else
cvx_begin quiet
end
variable x(n);
minimize( norm(A*x,2) );
subject to
norm([x;-sum(x)],1) <= 2*ub; % 约束条件 norm([x;-sum(x)],inf) <="ub/2;" c*x="=f;" cvx_end code>% 调用lcmv_l1per.m (in optimize_currents.m)
C = zeros(3*numOfROI,M);
f = zeros(3*numOfROI,1);
for n = 1:numOfROI
C(((n-1)*3+1):n*3,:) = [mean(A(tar_nodes{n},:),1);mean(A(tar_nodes{n}+Nlocs,:),1);mean(A(tar_nodes{n}+2*Nlocs,:),1);];
f(((n-1)*3+1):n*3) = [mean(d(tar_nodes{n}));mean(d(tar_nodes{n}+Nlocs));mean(d(tar_nodes{n}+2*Nlocs));];
end
status = '';
while strcmp(status,'Solved')~=1 %如果指定场强下无解，则降低场强
fprintf('CVX optimization infeasible under LCMV with L1 constraint on individual electrode...reducing target intensity...\n')
[s_opt,status] = lcmv_l1per(S*V',C,f,S_max,verbose);
f = f/1.25;
end
if all(abs(f)<1e-5 % 降到很低也没有结果 s_opt="nan(M,1);" error('cvx optimization infeasible for this target using lcmv with l1 constraint on individual electrode! consider targeting a nearby location.') end x_opt="A*s_opt;" < code>Optimizing for intensity （MI）
聚焦的优化对应于最大限度地提高tDCS的安全性，即不期望的大脑区域不会受到刺激。在某些情况下，可能需要牺牲聚焦性，而是最大化目标位置的电场。（4）的框架可以很容易地描述优化强度的问题：从（13）中回想，$C$由对应于目标节点的$A$的行组成，而$e_o$表示目标处的期望场方向。因此，目标处所需方向上的强度最大化采用线性规划问题的形式：
$e_0^TCs$是目标节点处的电场位于指向所需电场方向的矢量上。
程序中限定了默认最大场强法的电极数为4
% 最大强度法程序实现（max_l1per.m）
function [x,cvx_status] = max_l1per(f,ub,solElecNum,verbose)
n = size(f,1);
if verbose
cvx_begin
else
cvx_begin quiet
end
variable x(n);
maximize( sum(f'*x) );
subject to
norm([x;-sum(x)],1) <= 2*ub; norm([x;-sum(x)],inf) <="2*ub/solElecNum;" cvx_end code>% 调用max_l1per.m(in optimize_currents.m)
Cf = zeros(M,numOfROI);
for n = 1:numOfROI
Cx = mean(A(tar_nodes{n},:),1);
Cy = mean(A(tar_nodes{n}+Nlocs,:),1);
Cz = mean(A(tar_nodes{n}+2*Nlocs,:),1); % 将目标区域的A矩阵值在x,y,z上平均
C = [Cx;Cy;Cz];
f = [mean(d(tar_nodes{n}));mean(d(tar_nodes{n}+Nlocs));mean(d(tar_nodes{n}+2*Nlocs))]; % 将目标区域的期望场强在x,y,z做平均
Cf(:,n) = C'*f; % 最大化平均值
end
[s_opt,status] = max_l1per(Cf,S_max,solElecNum,verbose);
x_opt = A*s_opt;
