método multipolar rápido de elementos de contorno basado en carga

La formulación basada en carga del método de elementos de contorno (BEM) es una técnica numérica de reducción de dimensionalidad que se utiliza para modelar fenómenos electromagnéticos cuasiestáticos en medios conductores altamente complejos (como el cerebro humano ) con un número muy grande (hasta aproximadamente mil millones) de incógnitas. El BEM basado en carga resuelve una ecuación integral de la teoría del potencial ^{[ 1 ]} escrita en términos de la densidad de carga superficial inducida . Esta formulación se combina naturalmente con la aceleración del método multipolar rápido (FMM), y el método completo se conoce como BEM-FMM basado en carga . La combinación de BEM y FMM es una técnica común en diferentes áreas del electromagnetismo computacional y, en el contexto del bioelectromagnetismo, proporciona mejoras con respecto al método de elementos finitos . ^{[ 2 ] [ 3 ] [ 4 ]}

Junto con el BEM basado en potencial eléctrico más común, ^{[ 5 ] [ 6 ]} el BEM basado en carga cuasiestática, derivado en términos de la densidad de monocapa (carga), para un medio de un solo compartimento se conoce en la teoría del potencial ^{[ 1 ]} desde principios del siglo XX. Para medios conductores de múltiples compartimentos, la formulación de la densidad de carga superficial apareció por primera vez en forma discretizada (para interfaces facetadas) en el artículo de 1964 de Gelernter y Swihart. ^{[ 7 ]} Una forma continua posterior, que incluye efectos dependientes del tiempo y dieléctricos, apareció en el artículo de 1967 de Barnard, Duck y Lynn. ^{[ 8 ]} El BEM basado en carga también se ha formulado para medios conductores, dieléctricos y magnéticos, ^{[ 9 ]} y se ha utilizado en diferentes aplicaciones. ^{[ 10 ]}

En 2009, Greengard et al. ^{[ 11 ]} aplicaron con éxito el método BEM basado en carga con aceleración multipolar rápida a la electrostática molecular de dieléctricos. Un enfoque similar para el modelado realista del cerebro humano con múltiples compartimentos conductores fue descrito por primera vez por Makarov et al. ^{[ 12 ]} en 2018. Además, el método multipolar rápido multinivel basado en BEM se ha utilizado ampliamente en estudios de radar y antenas en frecuencias de microondas ^{[ 13 ]} , así como en acústica. ^{[ 14 ] [ 15 ]}

El método BEM basado en carga se basa en el concepto de un campo eléctrico impreso (o primario).${\displaystyle \mathbf {E} ^{i}}$y un campo eléctrico secundario${\displaystyle \mathbf {E} ^{s}}$El campo aplicado suele conocerse de antemano o es trivial de encontrar. Para el cerebro humano, el campo eléctrico aplicado se puede clasificar de la siguiente manera:

1. Un campo conservador${\displaystyle \mathbf {E} ^{i}}$derivado de una densidad impresa de fuentes de corriente EEG o MEG en un medio infinito homogéneo con la conductividad${\displaystyle \sigma }$en la ubicación de origen; ^{[ 16 ]}
2. Un campo solenoidal instantáneo${\displaystyle \mathbf {E} ^{i}}$de una bobina de inducción obtenida a partir de la ley de inducción de Faraday en un medio homogéneo infinito (aire), cuando se trata de problemas de estimulación magnética transcraneal (EMT); ^{[ 12 ] [ 17 ]}
3. Un campo superficial${\displaystyle \mathbf {E} ^{i}}$derivada de una densidad de corriente superficial impresa${\displaystyle \mathbf {J} ^{i}=\sigma \mathbf {E} ^{i}}$de electrodos de corriente que inyectan corriente eléctrica en el límite de un compartimento con conductividad${\displaystyle \sigma }$cuando se trata de estimulación transcraneal de corriente continua (tDCS) o estimulación cerebral profunda (DBS); ^{[ 18 ]}
4. Un campo conservador${\displaystyle \mathbf {E} ^{i}}$de cargas depositadas en electrodos de voltaje para tDCS o DBS. Este problema específico requiere un tratamiento acoplado ya que estas cargas dependerán del entorno; ^{[ 18 ]}
5. En aplicación al modelado multiescala , un campo${\displaystyle \mathbf {E} ^{i}}$obtenido a partir de cualquier otra solución numérica macroscópica en un dominio espacial pequeño (mesoescala o microescala) dentro del cerebro. Por ejemplo, se puede utilizar un campo constante. ^{[ 19 ]}

   

Cuando se "activa" el campo aplicado, las cargas libres ubicadas dentro de un volumen conductor D comienzan inmediatamente a redistribuirse y acumularse en los límites (interfaces) de regiones de diferente conductividad en D. Una densidad de carga superficial${\displaystyle \rho (\mathbf {r} )}$aparece en las interfaces de conductividad. Esta densidad de carga induce un campo eléctrico conservativo secundario.${\displaystyle \mathbf {E} ^{s}}$siguiendo la ley de Coulomb .

Un ejemplo es un ser humano bajo una línea eléctrica de corriente continua con un campo conocido. ${\displaystyle \mathbf {E} ^{i}}$Dirigida hacia abajo. La superficie superior del cuerpo conductor humano se cargará negativamente, mientras que su parte inferior se cargará positivamente. Estas cargas superficiales crean un campo eléctrico secundario que cancela o bloquea eficazmente el campo primario en todo el cuerpo, de modo que no fluye corriente dentro del cuerpo en condiciones de corriente continua constante.

Otro ejemplo es una cabeza humana con electrodos conectados. En cualquier interfaz de conductividad con un vector normal${\displaystyle \mathbf {n} }$ apuntando desde un compartimento "interior" (-) de conductividad${\displaystyle \sigma ^{-}}$a un compartimento "exterior" (+) de conductividad${\displaystyle \sigma ^{+}}$La ley de corrientes de Kirchhoff requiere la continuidad de la componente normal de la densidad de corriente eléctrica. Esto conduce a la condición de contorno interfacial en la forma

para cada faceta en una interfaz triangulada. Siempre que${\displaystyle \sigma ^{\pm }}$son diferentes entre sí, los dos componentes normales del campo eléctrico,${\displaystyle \mathbf {E} ^{\pm }\cdot \mathbf {n} }$También debe ser diferente. Tal salto a través de la interfaz solo es posible cuando existe una capa de carga superficial en dicha interfaz. Por lo tanto, si se aplica una corriente o voltaje eléctrico, la densidad de carga superficial se ve afectada.

El objetivo del análisis numérico es encontrar la distribución desconocida de carga superficial y, por lo tanto, el campo eléctrico total.${\displaystyle \mathbf {E} =\mathbf {E} ^{i}+\mathbf {E} ^{s}}$(y el potencial eléctrico  total si fuera necesario) en cualquier lugar del espacio.

A continuación, se presenta una derivación basada en la ley de Gauss y la ley de Coulomb . Todas las interfaces de conductividad, denotadas por S , se discretizan en facetas triangulares planas.${\displaystyle t_{m}}$ con centros${\displaystyle \mathbf {r} _{m}}$. Supongamos que una faceta m -ésima con el vector normal${\displaystyle \mathbf {n} _{m}}$ y área${\displaystyle A_{m}}$posee una densidad de carga superficial uniforme${\displaystyle \rho _{m}}$Si existiera una malla tetraédrica volumétrica, las facetas cargadas pertenecerían a tetraedros con diferentes valores de conductividad. Primero calculamos el campo eléctrico.${\displaystyle \mathbf {E} _ {m}^{+}}$en ese punto${\displaystyle \mathbf {r} _{m}+\delta \mathbf {n} _{m}}$, para${\displaystyle \delta \rightarrow 0^{+}}$es decir, justo fuera de la faceta 𝑚 en su centro. Este campo contiene tres contribuciones:

• El campo eléctrico impreso continuo${\displaystyle \mathbf {E} ^{i}}$sí mismo;
• Un campo eléctrico de la m -ésima faceta cargada misma. Muy cerca de la faceta, se puede aproximar como el campo eléctrico de una lámina infinita de carga superficial uniforme.${\displaystyle \rho _{m}}$. ^{[ 20 ]} Según la ley de Gauss , está dado por${\displaystyle +\rho _{m}/2\varepsilon _{0}\cdot \mathbf {n} _{m}}$dónde${\displaystyle \varepsilon _{0}}$es una permitividad eléctrica de fondo;
• Un campo eléctrico generado por todas las demás facetas${\displaystyle t_{n}}$, que aproximamos como cargas puntuales de carga${\displaystyle A_{n}\rho _{n}}$en cada centro${\displaystyle \mathbf {r} _{n}}$.

Un tratamiento similar se aplica al campo eléctrico.${\displaystyle \mathbf {E} _ {m}^{-}}$ justo dentro de la faceta 𝑚, pero el campo eléctrico de la lámina plana de carga cambia de signo. Usando la ley de Coulomb para calcular la contribución de facetas diferentes de${\displaystyle t_{m}}$, encontramos

De esta ecuación, vemos que la componente normal del campo eléctrico efectivamente experimenta un salto a través de la interfaz cargada. Esto es equivalente a una relación de salto de la teoría del potencial. ^{[ 1 ]} Como segundo paso, las dos expresiones para${\displaystyle \mathbf {E} _ {m}^{\pm }}$se sustituyen en la condición de contorno interfacial${\displaystyle \sigma ^{-}\mathbf {E} _{m}^{-}\cdot \mathbf {n} _{m}=\sigma ^{+}\mathbf {E} _{m}^{+}\cdot \mathbf {n} _{m}}$, aplicado a cada faceta 𝑚. Esta operación conduce a un sistema de ecuaciones lineales para densidades de carga desconocidas.${\displaystyle \rho _{m}}$lo cual resuelve el problema:

dónde${\displaystyle K_{m}={\frac {\sigma ^{-}-\sigma ^{+}}{\sigma ^{-}+\sigma ^{+}}}}$es el contraste de conductividad eléctrica en la faceta m . La constante de normalización${\displaystyle \varepsilon _{0}}$se cancelará después de que la solución se sustituya en la expresión para${\displaystyle \mathbf {E} ^{s}}$y se vuelve redundante.

Para caracterizaciones modernas de topologías cerebrales con niveles de complejidad cada vez mayores, el sistema de ecuaciones anterior para${\displaystyle \rho _{m}}$es muy grande; por lo tanto, se resuelve iterativamente. Una estimación inicial para${\displaystyle \rho _{m}}$es el último término en su lado derecho mientras que la suma se ignora. A continuación, se calcula la suma y se refina la estimación inicial, etc. Esta solución ^{[ 12 ] [ 21 ]} emplea el método iterativo de Jacobi simple. El método generalizado de residuos mínimos (GMRES) , más riguroso, produce una convergencia mucho más rápida del BEM-FMM. ^{[ 2 ] [ 3 ] [ 16 ] [ 17 ] [ 18 ]} En cualquier caso, el trabajo principal consiste en calcular la suma entre corchetes en el sistema de ecuaciones anterior para cada${\displaystyle {m}}$en cada iteración; esta operación corresponde a una multiplicación repetitiva matriz-vector. Sin embargo, se puede reconocer esta suma como un campo eléctrico (veces${\displaystyle {\frac {1}{2\pi \epsilon _ {0}}}}$) de${\displaystyle {M}}$cargos a calcular en${\displaystyle {M}}$ puntos de observación. Tal cálculo es precisamente la tarea del método multipolar rápido , que realiza una multiplicación rápida de matriz por vector en${\displaystyle O(M\log {M})}$ o incluso${\displaystyle O(M)}$ operaciones en lugar de${\displaystyle O(M^{2})}$La biblioteca FMM3D ^{[ 22 ]} , implementada tanto en Python como en MATLAB , puede utilizarse para este propósito. Por lo tanto, no es necesario formar ni almacenar la matriz densa del sistema típica del BEM estándar .

El sistema de ecuaciones formulado anteriormente se deriva mediante el método de colocación y es menos preciso. ^{[ 11 ]} La ecuación integral correspondiente se obtiene a partir de las relaciones de salto locales de la teoría del potencial ^{[ 23 ]} y la condición de contorno interfacial local de continuidad de la corriente eléctrica normal. Es una ecuación integral de Fredholm de segundo tipo.

Su derivación no involucra las identidades de Green (integraciones por partes) y es aplicable a geometrías no anidadas. Cuando se aplica el método de Galerkin y se siguen utilizando las mismas funciones base de orden cero (con una densidad de carga constante para cada faceta) en interfaces trianguladas, obtenemos exactamente la misma discretización que antes si reemplazamos las integrales dobles sobre superficies.${\displaystyle S_{m}}$y${\displaystyle S_{n}}$de triángulos${\displaystyle t_{m}}$y${\displaystyle t_{n}}$, respectivamente, por

${\displaystyle \int _{S_{m}}\int _{S_{n}}{\frac {\mathbf {r} -\mathbf {r^{\prime }} }{|\mathbf {r} -\mathbf {r^{\prime }} |^{3}}}ds(\mathbf {r^{\prime }} )ds(\mathbf {r} )\approx {A_{m}}{A_{n}}{\frac {\mathbf {r} _{m}-\mathbf {r} _{n}}{|\mathbf {r} _{m}-\mathbf {r} _{n}|^{3}}},}$

dónde${\displaystyle {A_{n}}}$es el área de la superficie del triángulo${\displaystyle {t_{n}}}$Esta aproximación solo es válida cuando${\displaystyle |\mathbf {r} _{m}-\mathbf {r} _{n}|}$es mucho mayor que el tamaño típico de una faceta, es decir, en el "campo lejano". De lo contrario, se deben usar  fórmulas semianalíticas ^{[ 24 ] [ 25 ]} y cuadraturas gaussianas para triángulos ^{[ 26 ] [ 12 ]} . Normalmente, se deben precalcular, almacenar y usar de 4 a 32 integrales vecinas por faceta en cada iteración. ^{[ 12 ] [ 2 ] [ 17 ] [ 18 ] [ 27 ]} Esta es una corrección importante al método multipolar rápido simple en el "campo cercano" que también debe usarse en la formulación discreta simple derivada anteriormente. Dicha corrección permite obtener una resolución numérica (pero no anatómica) sin restricciones en el cerebro. ^{[ 17 ]}

Las aplicaciones del BEM-FMM basado en carga incluyen el modelado de la estimulación cerebral ^{[ 3 ] [ 17 ] [ 18 ] [ 21 ]} con cálculos de TMS precisos casi en tiempo real ^{[ 28 ] [ 4 ]} , así como registros neurofisiológicos. ^{[ 16 ]} También incluyen el modelado de topologías de cabeza mesoescala desafiantes, como membranas cerebrales delgadas ^{[ 29 ] [ 27 ]} ( duramadre , aracnoides y piamadre ). Esto es particularmente importante para predicciones precisas de la dosis de estimulación transcraneal de corriente continua y terapia electroconvulsiva . ^{[ 30 ]} El BEM-FMM permite un refinamiento de malla adaptativo sencillo que incluye múltiples compartimentos cerebrales extracerebrales. ^{[ 27 ] [ 29 ]} Otra aplicación es el modelado de perturbaciones del campo eléctrico dentro de un árbol neuronal/axonal densamente empaquetado. ^{[ 19 ]} Dichas perturbaciones cambian la función de activación biofísica . Se está desarrollando una formulación BEM basada en carga para un prometedor modelado biofísico bidominio de procesos axonales. ^{[ 31 ]}

En su forma actual, el BEM-FMM basado en carga solo es aplicable a medios homogéneos por partes multicompartimentales; no puede manejar tejidos macroscópicamente anisotrópicos . Además, el número máximo de facetas (grados de libertad) está limitado a aproximadamente${\displaystyle 10^{9}}$para los recursos de hardware informático académicos típicos utilizados a partir de 2023.

• Electromagnetismo computacional
• Método de elementos de contorno
• Método multipolar rápido
• Neurociencia computacional
• Estimulación magnética transcraneal
• Estimulación transcraneal de corriente continua
• Electroencefalografía
• Magnetoencefalografía

• Un estudio sobre ecuaciones integrales para el modelado bioeléctrico , preimpresión .
• Sitio del proyecto FMM3D de la Fundación Simons en GitHub del Instituto Flatiron .