# A mathematical model for coregistered data from electroencephalography and diffusive optical tomography

M. Darbas, S. Lohrengel, B. Sulis

► To cite this version:

M. Darbas, S. Lohrengel, B. Sulis. A mathematical model for coregistered data from electroencephalography and diffusive optical tomography. Mathematical Modelling of Natural Phenomena, 2025, 20, pp.4. 10.1051/mmnp/2025001. hal-04943895

HAL Id: hal-04943895

https://hal.science/hal-04943895v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A MATHEMATICAL MODEL FOR COREGISTERED DATA
FROM ELECTROENCEPHALOGRAPHY AND DIFFUSIVE
OPTICAL TOMOGRAPHY

M. DARBAS<fnref n="1" />, S. LOHRENGEL<fnref n="2" />,* AND B. SULIS<fnref n="2" />

**Abstract.** A mathematical model for coregistered time-dependent electroencephalography (EEG) and diffusive optical tomography (DOT) is developed and analysed. Evolution with time is introduced by considering time-dependent dipolar sources in the EEG model and time-dependent optical parameters for DOT. Dimensional analysis shows that time-derivatives can be neglected. A non-linear system of differential equations from literature is used to model the postsynaptic current and hemodynamic parameters at the neuron level. A key point of the full model is to explain how these quantities provide, at the level of the whole head, the moment of the dipolar source term of the EEG problem and the behaviour in time of the optical parameters of the DOT model. The well-posedness of the time-dependent EEG problem is proved by the subtraction approach for moments with $L^2$-regularity in time and continuous source trajectories. For the time-dependent DOT model with continuous optical parameters in time, standard results of functional analysis apply. We explain the full pipeline from the stimulation current up to the simulated signals recorded at the electroptodes. Numerical results for a three-dimensional realistic head model illustrate the capacity of simultaneous EEG/DOT measurements to attest neurovascular coupling between the neural activity and changes in the hemodynamic parameters.

**Mathematics Subject Classification.** 92C50, 92-10, 65N30.

Received May 28, 2024. Accepted January 9, 2025.

# 1. INTRODUCTION

The challenges of understanding the human brain are major and of great interest to better diagnosis and cure of many brain diseases. Various brain imaging modalities have been developed over the years, each dedicated to the investigation of some specific brain functions. We can cite functional magnetic resonance imaging (fMRI), computerized tomography (CT), positron emission tomography (PET), electroencephalography (EEG), magnetencephalography (MEG), functional near-infrared spectroscopy (fNIRS) or diffuse optical tomography (DOT). Some of them provide complementary information on the brain and its disorders, and the benefits of multimodal neuroimaging are studied in the last decades. In this paper, we focus on the synchronous measurements of EEG and DOT.

<footnote>*Keywords and phrases:* Time-dependent electroencephalography with dipolar sources, time-dependent diffusive optical tomography, neurovascular coupling, neuronal model.</footnote>

<footnote n="1">¹ LAGA UMR CNRS 7539, Université Sorbonne Paris Nord, France.</footnote>

<footnote n="2">² LMR UMR CNRS 9008, Université de Reims-Champagne Ardenne, Reims, France.</footnote>

<footnote>* Corresponding author; stephanie.lohrengel@univ-reims.fr</footnote>
On the one hand, electroencephalography (EEG) is one of the most widespread functional brain imaging techniques. Measurements of the electric potential generated by normal or pathological brain activity are taken at electrodes attached to the scalp. They record in a passive and non-invasive way the voltage potential fluctuations between different cortical regions. EEG-monitoring can be done at the bedside of the patient which makes the technique particularly appropriate for neonates and premature babies. The important goal of brain imaging using EEG is to localize cerebral sources generating measured EEG signals. EEG is known to have an excellent temporal resolution and is able to record neural events in order of one millisecond. Its spatial resolution is generally limited (a few centimeters). Neonatal EEG is used to assess seizure recognition and classification and to make epilepsy syndrome diagnosis (e.g. [1]). It provides also prognostic information for other brain dysfunctions (e.g. [2]).

On the other hand, diffuse optical tomography (DOT) is a functional brain imaging technique that measures at optodes (sensors) on the scalp the hemodynamic changes in the brain. DOT is an extension of fNIRS (see e.g. [3]). More precisely, DOT aims in localizing changes in oxy- and deoxyhemoglobin (that result from brain activity) within the brain by using absorption of near-infrared light. Similar to EEG, it is a non-invasive technique. DOT may be seen as an alternative to fMRI. In comparison to fMRI, the DOT system has the advantage to be relatively low-cost, portable and bedside-compatible.

The development of multimodal analysis yields a more detailed understanding of brain dynamics. We are particularly interested in the coupling of EEG with optical imaging, taking advantage of the possible coregistration of EEG and NIRS which has already been operated by collaborators from GRAMFC INSERM UMR-S 1105 (Amiens' hospital) [4] for neonates and premature babies. EEG measures directly neuronal activity while optical imaging techniques (DOT, NIRS, fNIRS) capture the physiological changes linked to neurovascular coupling. To sum up, in response to a certain task (due to normal or pathological origin), the neural activity increases in a brain area. A metabolic demand, *i.e.* changes in oxy- and deoxyhemoglobin, is observed generally in that area. Coregistered data from EEG and DOT give information to clinicians for understanding the neurovascular coupling which occurs in some brain disorders. The coupling of the two modalities is more particularly used in pediatric clinical research [5, 6]. An example is the study of [7] where hemodynamic changes have been observed among 10s before an epileptic seizure. Another clinical investigation using simultaneous EEG and DOT is described in [8].

The mathematical modelling of EEG and DOT has been widely addressed in the literature. Without being exhaustive, we cite [9–12] and references therein for an analysis of the forward EEG problem where the models are static. Experimental DOT devices mainly rely on three different types of light sources: systems with constant light source measure the total intensity of light at the surface of the object, systems with frequency modulated sources operate at radio-frequency and allow to measure separately the amplitude as well as the phase shift of the intensity, and time resolved systems are based on sources of an ultra short time impulse. Time-resolved DOT is an existing imaging modality [13] that is actually less relevant for clinical purposes due to the high cost of the technical device. In this paper, we thus limit our investigations for DOT with frequency modulated source terms. We refer to [14–16] for an overview of models for DOT based on the diffusion approximation which are time-harmonic in the case of frequency modulated sources with stationary optical parameters.

The present work is dedicated to the development and the analysis of a mathematical model for simulating coregistered time-dependent EEG-DOT data. Given an input stimulus, the underlying forward problem consists in computing the electric potential and the photon density which are generated by the brain activity corresponding to the stimulus, at electroptodes on the scalp. To this end, we combine different models at the neuron level and at the level of the whole head. First, we need to model the neurovascular coupling mechanism, that is the dynamics of the neural activity and its vascular response. This is done by solving a system of non-linear ordinary differential equations. Different models have been developed in the past (see e.g. [17–22] and references therein). In this paper, we focus on a model for neuronal activity based on the work of [21] coupled to the balloon model [18, 19] for the hemodynamic response. These models are used to get the behaviour in time of the postsynaptic current and the hemodynamic parameters at the neuron level. Next, we search space-time models that allow to propagate and diffuse this local brain information at the macroscopic level of the whole head. Full Maxwell's equations and the Radiative Transfer Equation (RTE) are computationally
expensive. Instead, we introduce time evolution in the existing stationary PDE models for EEG and DOT by means of time-dependent source terms for EEG and time-dependent optical parameters for DOT. The absence of time-derivatives is justified by dimensional analysis. The behaviour in time of the EEG source term and the optical coefficients for DOT are deduced from the ODE's models at the neuron level. Finally, we get a complete pipeline of models that allow the numerical simulation of coupled EEG-DOT signals from a (time-dependent) input stimulus up to the electric potential and photon density which depend on space and time variables.

The paper is organized as follows. In Section 2, we justify a time-dependent model for the EEG forward problem and we address an existence and uniqueness result using the subtraction method. In Section 3, we propose a model for the diffusion approximation with time-dependent optical parameters in a frequency-modulated setting. Section 4 is devoted to the modelling of the neurovascular coupling at the neuron level. Finally, we describe the whole numerical pipeline in Section 5. Simulations of coregistered EEG/DOT signals for a realistic three-dimensional head model are performed.

## 2. A MATHEMATICAL MODEL FOR TIME-DEPENDENT EEG

### 2.1. Dimensional analysis of Maxwell's equations

We study an electromagnetic phenomenon in a space-time domain with characteristic length $\ell$ and time scale $\tau$. The continuous medium is characterized by its electric permittivity $\varepsilon$, magnetic permeability $\mu$ and conductivity $\sigma$. Maxwell's equations describe the propagation of an electromagnetic field in this medium. Under the assumption that the constitution laws of the medium are linear, they are given by

$$ 
\begin{aligned}
\nabla \cdot \mathbf{B} &= 0, \\
\nabla \times \mathbf{E} &= -\partial_t \mathbf{B},
\end{aligned} 
\qquad 
\begin{aligned}
\nabla \times \mathbf{B} &= \mu(\varepsilon\partial_t\mathbf{E} + \mathbf{J}), \\
\nabla \cdot \mathbf{E} &= \frac{\rho}{\varepsilon},
\end{aligned} 
\qquad (2.1) 
$$

and model the interaction between the electric field $\mathbf{E}$, the magnetic induction $\mathbf{B}$, the charge density $\rho$ and the current density $\mathbf{J}$. In a conducting medium with conductivity $\sigma$, Ohm's law states that the free current density is given by $\mathbf{J} = \sigma\mathbf{E}$. In the presence of impressed current sources, a source term $\mathbf{j}^p$ should be added to $\mathbf{J}$.

The static limit of Maxwell's equations is well understood in the case where all fields and sources are time-independent. As soon as there is some time-dependence, however, the full Maxwell system should apply unless mathematical analysis shows that some terms can be neglected. EEG measures the cerebral activity during the observation. This activity is clearly variable in time, meaning that time-dependent source terms should be considered. In this context, dimensional analysis is an interesting tool to reduce the set of Maxwell's equations and get simplified models under particular hypotheses. The idea is to observe the different physical quantities in terms of their units or dimensions. Smallness of some of them allows to neglect certain coupling between electric and magnetic fields [23]. To this end, we introduce the velocity of the system with modulus $v = \ell/\tau$. The light celerity in the medium with electromagnetic parameters $\varepsilon$ and $\mu$ is given by $c = 1/\sqrt{\varepsilon\mu}$. Here, we are interested in limit configurations where $v \ll c$.

For going further into details in the analysis, a first step is to express the quantities $\ell$ and $c$ in terms of the electromagnetic parameters $\varepsilon$, $\mu$ and $\sigma$. Following [23], we introduce the quantity $\tau_e = \varepsilon/\sigma$ which has dimension of time and can be interpreted as the electric charge diffusion time. Similarly, the quantity $\ell c\mu\sigma$ can be shown to be dimensionless and can be written as the quotient of $\tau_m = \mu\sigma\ell^2$ and $\tau_{em} = \ell/c$ which both have dimension of time. We next proceed to the scaling of Maxwell's equations. To this end, set $\mathbf{E}(t, \mathbf{x}) = e\mathbf{E}'(t', \mathbf{x}')$ and $\mathbf{B}(t, \mathbf{x}) = b\mathbf{B}'(t', \mathbf{x}')$ where $e, b$ are reference quantities and $\mathbf{E}', \mathbf{B}'$ are dimensionless quantities of order $\mathcal{O}(1)$ depending on $t' = t/\tau$ and $\mathbf{x}' = \mathbf{x}/\ell$. In the sequel, prime notation corresponds to dimensionless operators or variables.
Using the definition of the characteristic times, Faraday's and Ampère's laws $\nabla \times \mathbf{E} = -\partial_t \mathbf{B}$ and $\nabla \times \mathbf{B} = \mu(\varepsilon\partial_t\mathbf{E} + \sigma\mathbf{E})$ become respectively

$$ \nabla' \times \mathbf{E}' = -\frac{\tau_{em}}{\tau} \frac{cb}{e} \partial_{t'} \mathbf{B}', \qquad (2.2) $$

and

$$ \nabla' \times \mathbf{B}' = \frac{\tau_{em}}{\tau} \frac{e}{cb} \partial_{t'} \mathbf{E}' + \frac{\tau_{em}}{\tau_e} \frac{e}{cb} \mathbf{E}'. \qquad (2.3) $$

We are interested in the low frequency range where the characteristic dimensions of the system are such that

$$ \tau_{em} \ll \tau, \quad \tau_e \ll \tau, \quad \text{and} \quad \tau_m \ll \tau. \qquad (2.4) $$

Then, we cannot have simultaneously $\frac{\tau_{em}}{\tau} \frac{cb}{e} = O(1)$ and $\frac{\tau_{em}}{\tau} \frac{e}{cb} = O(1)$ and at least one of the time derivatives in (2.2) or (2.3) has to be neglected. We aim to determine the limit model from the only order between the characteristic time scales.

One may notice that $\tau_{em}^2 = \tau_e\tau_m$ which amounts to saying that any order between $\tau_{em}$ and $\tau_e$ induces an order of $\tau_m$ with respect to $\tau_{em}$. We distinguish three cases: $\tau_{em} \ll \tau_e$, $\tau_e \ll \tau_{em}$ and $\tau_e \sim \tau_{em}$. We choose to only detail the first one. According to assumption (2.4) on $\tau_e$, the relation $\tau_{em} \ll \tau_e$ implies the order $\tau_m \ll \tau_{em} \ll \tau_e \ll \tau$. From $\tau_e \ll \tau$, we deduce that $\frac{\tau_{em}}{\tau} \ll \frac{\tau_{em}}{\tau_e}$. The displacement current in the right hand side of (2.3) can thus be neglected and we get

$$ \nabla \times \mathbf{B} = \mu\mathbf{J}. \qquad (2.5) $$

Since the left hand side of (2.3) is of order 1, we further have $\frac{cb}{e} \sim \frac{\tau_{em}}{\tau_e} \ll 1$. Hence, the right hand side in (2.2) can be neglected, too, which yields

$$ \nabla \times \mathbf{E} = 0. \qquad (2.6) $$

The cases $\tau_e \ll \tau_{em}$ and $\tau_e \sim \tau_{em}$ lead to the same conclusions. Thus, in a frequency range satisfying (2.4), we get the following approximation of Maxwell's equations

$$ \begin{aligned} \nabla \cdot \mathbf{B} &= 0, & \nabla \times \mathbf{B} &= \mu\mathbf{J}, \\ \nabla \times \mathbf{E} &= 0, & \nabla \cdot \mathbf{E} &= \frac{\rho}{\varepsilon}, \end{aligned} \qquad (2.7) $$

called the quasi-stationary model.

## 2.2. Quasi-stationary approximation for EEG with time-dependent sources

In this section, we compare the order of magnitude of the characteristic times for the electromagnetic parameters of the different head tissues. In mathematical modelling of brain activity, one assumes that the magnetic permeability of head tissues is close to the one of free space. The parameters of the medium are given by $\varepsilon = \varepsilon_r\varepsilon_0$ and $\mu = \mu_0$ where $\varepsilon_0 = 8.85 \times 10^{-12} \text{F.m}^{-1}$ and $\mu_0 = 4\pi \times 10^{-7} \text{H.m}^{-1}$ are respective vacuum permittivity and permeability. The relative permittivity $\varepsilon_r$ and the conductivity $\sigma$ depend on the tissue, but also on the frequency range. The typical frequency range of the electromagnetic fields in the head is inferior to 100 Hz [24] and we thus take $\tau \approx 0.01$s as the characteristic time scale. We consider a multilayer head model which distinguishes
<table><caption>TABLE 1. Electric parameter set at 100 Hz [25, 28] and characteristic times.</caption><thead><tr><th></th><th>White matter</th><th>Grey matter</th><th>CSF</th><th>Skull</th><th>Scalp</th></tr></thead><tbody><tr><td>ε<sub>r</sub></td><td>1.67 × 10<sup>6</sup></td><td>3.91 × 10<sup>6</sup></td><td>1.1 × 10<sup>2</sup></td><td>5.85 × 10<sup>3</sup></td><td>1.14 × 10<sup>3</sup></td></tr><tr><td>σ (S m<sup>-1</sup>)</td><td>0.14</td><td>0.33</td><td>1.8</td><td>0.04</td><td>0.33</td></tr><tr><td>τ<sub>m</sub> (s)</td><td>2.53 × 10<sup>-9</sup></td><td>5.97 × 10<sup>-9</sup></td><td>3.25 × 10<sup>-8</sup></td><td>7.2 × 10<sup>-10</sup></td><td>6.0 × 10<sup>-9</sup></td></tr><tr><td>τ<sub>em</sub> (s)</td><td>5 × 10<sup>-7</sup></td><td>8 × 10<sup>-7</sup></td><td>4.19 × 10<sup>-9</sup></td><td>3.1 × 10<sup>-8</sup></td><td>1.35 × 10<sup>-8</sup></td></tr><tr><td>τ<sub>e</sub> (s)</td><td>1.05 × 10<sup>-5</sup></td><td>1.05 × 10<sup>-4</sup></td><td>5.4 × 10<sup>-10</sup></td><td>1.3 × 10<sup>-6</sup></td><td>3.1 × 10<sup>-8</sup></td></tr></tbody></table>

between white and grey matter and takes into account the cerebrospinal fluid (CSF) as well as skull and scalp. Uncertainty about the electric parameters should be taken into account since *in vivo* measurements in the human body are in general not available, especially for neonates and premature babies. In Table 1, we take the values for the electric permittivity from the IT'IS data base [25] which are given for a specific frequency from a parametrization fit based on the dispersion relation owing to [26]. The conductivity values are taken from [27, 28] for neonates. We recall the formulæ for the three characteristic times

$$ \tau_e = \varepsilon/\sigma, \tau_m = \mu\sigma\ell^2, \tau_{em} = \sqrt{\tau_e\tau_m}. $$

The characteristic length scale $\ell$ should be given by the thickness of the tissue in consideration. Notice, however, that $\tau_e$ is independent from the characteristic length $\ell$ whereas $\tau_m \sim \ell^2$ and $\tau_{em} \sim \ell$. Hence, if the quasi-static approximation is valid for a given $\ell$, it is also valid for any smaller length scale. We thus take for $\ell$ the dimensions of the neonatal head, i.e. $\ell \approx 0.12$ m. We then deduce the values of the three characteristic times.

From Table 1, we infer $\tau_m \ll \tau_{em} \ll \tau_e \ll \tau$ for white matter, grey matter and skull. In CSF, the order is reversed and in the scalp, the three characteristic times are approximately of the same order with $\tau_m \sim \tau_{em} \sim \tau_e \ll \tau$. These situations correspond to the different cases studied in the previous dimensional analysis. Thus, for modelling EEG, the quasi-stationary approximation (2.7) of the full Maxwell equations can be considered as a valid model. We focus on values for neonates and premature babies, but the conclusions of the dimensional analysis hold true for adults as well, the main difference being in the value of $\sigma$ in the skull which is generally assumed to be ten times smaller than in neonates.

Equation $\nabla \times \mathbf{E} = 0$ shows that the electric field $\mathbf{E}$ derives from a scalar electric potential, denoted by $u$, and is written as

$$ \mathbf{E}(t, \mathbf{x}) = -\nabla u(t, \mathbf{x}), \quad (t, \mathbf{x}) \in (0, T) \times \Omega, \qquad (2.8) $$

with $T > 0$ an observation time and $\Omega \subset \mathbb{R}^d$ ($d = 2, 3$) a head model. In the brain and in particular in the cortex, the synchronized effect among a multitude of neurons creates an intracellular current denoted by $\mathbf{j}^p$. The current density $\mathbf{J}$ produced by cerebral activity thus splits into two terms

$$ \mathbf{J} = \sigma \mathbf{E} + \mathbf{j}^p. \qquad (2.9) $$

By replacing (2.8) and (2.9) in the equation $\nabla \times \mathbf{B} = \mu \mathbf{J}$, and by applying the divergence operator, we obtain

$$ 0 = \nabla \cdot (\nabla \times \mathbf{B}) = \mu \nabla \cdot (\mathbf{j}^p - \sigma \nabla u), \qquad (2.10) $$

which gives the following elliptic equation

$$ \nabla \cdot (\sigma \nabla u) = \nabla \cdot \mathbf{j}^p. \qquad (2.11) $$
In order to numerically reproduce spatially localized phenomena which are variable in time, we propose a time-dependent model for the source term $\mathbf{j}^p$ in the form of a set of $M$ electrical current dipoles located in the brain

$$ \mathbf{j}^p(t, \cdot) = \sum_{m=1}^{M} \mathbf{q}_m(t)\delta(\cdot - S_m(t)), \qquad (2.12) $$

where $\delta(\cdot)$ denotes the Dirac delta distribution. Here, $S_m(t) \in \Omega$ and $\mathbf{q}_m(t) \in \mathbb{R}^d$ are, respectively, the position and the moment of the $m$-th current dipole at time $t$ which is situated in the subdomain of $\Omega$ that models the brain or, more precisely, the grey matter. The vector $\mathbf{q}_m(t)$ defines the strength and direction of the $m$-th dipole at time $t$ and is expressed in (A m). Multiplication by the Dirac delta yields (A m$^{-2}$) for the unit of the current density $\mathbf{j}_p$. The right hand side of (2.11) then reads

$$ F(t, \cdot) := \nabla \cdot \mathbf{j}^p(t, \cdot) = \sum_{m=1}^{M} \mathbf{q}_m(t) \cdot \nabla \delta(\cdot - S_m(t)). \qquad (2.13) $$

Assuming that the conductivity of air is zero and that no electric current can flow out of the scalp, the electric potential $u$ is thus solution of the following boundary problem with homogeneous Neumann condition

$$ \left\{ \begin{array}{ll} \nabla \cdot (\sigma \nabla u) = F & \text{in } (0, T) \times \Omega, \\ \sigma \partial_n u = 0 & \text{on } (0, T) \times \partial \Omega, \end{array} \right. \qquad (2.14) $$

where $u$ and $F$ are functions of space and time. Problem (2.14) is the instationary extension of the standard forward problem for EEG.

## 2.3. The subtraction approach for EEG with time-dependent sources

In this section, we address the resolution of the forward problem (2.14). Mathematically, a head model can be described as follows. Let $\Omega \subset \mathbb{R}^d$ be a bounded simply connected domain with regular boundary $\Gamma := \partial \Omega$ and consider a partition of $\Omega$ into $P$ open subdomains $(\Omega_p)_{p=1,...,P}$, such that

$$ \bar{\Omega} = \bigcup_{p=1}^{P} \bar{\Omega}_p \text{ and } \Omega_p \cap \Omega_q = \emptyset \quad \forall p \neq q. $$

Subdomains $\Omega_p$ describe the different tissues of the head. In the case of concentric subdomains as in Figure 1, we denote by $\Gamma_p$ the interface between the subdomains $\Omega_p$ and $\Omega_{p+1}$ and assume that $(\Gamma_p)_p$ are closed regular surfaces. Let $\mathbf{n}_p$ be the unit normal vector to $\Gamma_p$ from $\Omega_p$ to $\Omega_{p+1}$.

Notice that (2.14) includes the following transmission conditions at any interface $\Gamma_p$ between two subdomains $\Omega_p$ and $\Omega_q$

$$ [u]|_{\Gamma_p} = [\sigma \partial_n u]|_{\Gamma_p} = 0 \text{ on } (0, T) \times \Gamma_p. \qquad (2.15) $$

Here, $[f]|_{\Gamma_p} = f|_{\Omega_q} - f|_{\Omega_p}$ denotes the jump across the interface $\Gamma_p$ of the function $f$ defined on $\Omega$.

We make the following assumptions on the moments, sources, and the conductivity:

(H1) $\mathbf{q}_m \in L^2(0, T)^d$ and $S_m \in C^0([0, T])^d, \forall m \in \{1, ..., M\}$.

(H2) At time $t$, the points $(S_m(t))_m$ are mutually distinct, i.e. $S_m(t) \neq S_k(t), \forall m \neq k$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Three-layer head model.</figcaption></figure>

(H_3) All sources are located in the same subdomain $\Omega_{p_0}$, $p_0 \in \{1, \dots, P\}$, and there is a domain $D \subset \subset \Omega_{p_0}$, such that

$$ \bigcup_{m=1}^{M} \{S_m(t), t \in [0, T]\} \subset D, \qquad (2.16) $$

and the conductivity is constant on $D$: $\sigma|_{D} = \sigma_{D}$ for a constant $\sigma_{D} > 0$.

Here, the choice of moments $q_m$ with $L^2$-regularity in time is motivated by a possible need of signal processing tools in the analysis of EEG measurements. From a mathematical point of view, signals are naturally assumed to belong to some subspace of $L^2$.

Due to the lack of regularity in space of the source term $F$, a direct variational formulation of (2.14) in $H^1(\Omega)$ is not possible. A possibility is to apply the subtraction approach [12, 29, 30]. It consists of decomposing the potential $u$ into a potential $\tilde{u}$ which contains the singularity and a regular lifting $w$:

$$ u = \tilde{u} + w \text{ on } (0, T) \times \Omega. $$

But in the time-dependent setting with moments $q_m$ belonging to $L^2(0, T)^d$, the boundary value problem (2.14) with source term (2.13) can not be understood at fixed time $t$. The definition of the singular potential $\tilde{u}$ thus needs to be done with caution. To this end, we introduce for $m \in \{1, \dots, M\}$ and $i \in \{1, \dots, d\}$, the canonical source term $F_m^{(i)}$:

$$ \forall t \in [0, T], \quad F_m^{(i)}(t, \cdot) = \mathbf{e}^{(i)} \cdot \nabla \delta(\cdot - S_m(t)) \qquad (2.17) $$

where $\mathbf{e}^{(i)}$ denotes the $i$-th canonical basis vector of $\mathbb{R}^d$. Notice that $F_m^{(i)}$ is well defined at any time step $t$ since $S_m$ is continuous on $[0, T]$. The associated singular potential $\tilde{u}_m^{(i)}(t, \cdot)$, $i = 1, \dots, d$, is solution of the following Poisson equation

$$ \forall t \in [0, T], \quad \sigma_D \Delta \tilde{u}_m^{(i)}(t, \cdot) = F_m^{(i)}(t, \cdot) \text{ in } \mathbb{R}^d. $$

Thus, $\tilde{u}_m^{(i)}(t, \cdot)$ is obtained by convolution in the space variable $\mathbf{x}$ of the fundamental solution of the Laplace equation with the right hand side $\frac{1}{\sigma_D} \mathbf{e}^{(i)} \cdot \nabla \delta(\cdot - S_m(t))$ which leads to

$$ \tilde{u}_m^{(i)}(t, \mathbf{x}) = \frac{1}{2^{d-1}\pi\sigma_D} \mathbf{e}^{(i)} \cdot \frac{\mathbf{x} - S_m(t)}{|\mathbf{x} - S_m(t)|^d}, \quad \forall t \in [0, T], \forall \mathbf{x} \in \mathbb{R}^d \setminus \{S_m(t)\}. \qquad (2.18) $$
At fixed time $t \in [0, T]$, consider the following boundary value problem with regular right hand side

$$ egin{cases} -\nabla \cdot (\sigma \nabla w_m^{(i)}(t, \cdot)) = \nabla \cdot ((\sigma - \sigma_D)\nabla \tilde{u}_m^{(i)}(t, \cdot)) & \text{in } \Omega, \\ \sigma \partial_n w_m^{(i)}(t, \cdot) = -\sigma \partial_n \tilde{u}_m^{(i)}(t, \cdot) & \text{on } \Gamma. \end{cases} \quad (2.19) $$

We prove here below that (2.19) admits a variational solution in an appropriated vector space. To this end, define the bilinear form $a(\cdot, \cdot)$ on $H^1(\Omega) \times H^1(\Omega)$ by

$$ a(w, v) = \int_{\Omega} \sigma \nabla w \cdot \nabla v \, dx. $$

For fixed $t \in [0, T]$, denote by $l_m^{(i)}(t; \cdot)$ the following linear form defined for $v \in H^1(\Omega)$,

$$ l_m^{(i)}(t; v) = \int_{\Omega \setminus \bar{D}} (\sigma_D - \sigma) \nabla \tilde{u}_m^{(i)}(t, \cdot) \cdot \nabla v \, dx - \int_{\Gamma} \sigma_D \partial_n \tilde{u}_m^{(i)}(t, \cdot) v \, ds. $$

Since problem (2.19) involves a Neumann boundary condition, its solution is determined up to an additive constant only. We therefore introduce the subspace of $H^1(\Omega)$ of functions with vanishing mean value,

$$ V = H^1(\Omega) \cap L_0^2(\Omega), \quad (2.20) $$

where

$$ L_0^2(\Omega) = \left\{ v \in L^2(\Omega) \mid \int_{\Omega} v \, dx = 0 \right\} $$

and $H^1(\Omega)$ denotes the standard Sobolev space of square integrable real-valued functions with derivatives up to the first order.

Then, the following proposition holds true:

**Proposition 2.1.** Let $\sigma \in L^\infty(\Omega)$ be such that $0 < \sigma_{\min} \le \sigma \le \sigma_{\max}$ a.e. on $\Omega$ with constants $\sigma_{\min}$ and $\sigma_{\max}$. Let $S_m \in C^0(0, T)^d$ for any $m \in \{1, \dots, M\}$ and assume that $(H_2)$ and $(H_3)$ are satisfied. Then, for $m \in \{1, \dots, M\}$ and $i \in \{1, \dots, d\}$, the following problem admits a unique solution:

$$ egin{cases} \text{Find } w_m^{(i)} \in C^0([0, T]; V) \text{ such that} \\ a(w_m^{(i)}(t, \cdot), v) = l_m^{(i)}(t; v) \quad \forall v \in H^1(\Omega), \forall t \in [0, T]. \end{cases} \quad (2.21) $$

The solution of (2.21) satisfies (2.19) in a weak sense.

*Proof.* According to the assumptions on $\sigma$, the bilinear form $a(\cdot, \cdot)$ is continuous and coercive on $V \times V$. For the linear form, consider $v \in H^1(\Omega)$ and $t \in [0, T]$. Since the potential $\tilde{u}_m^{(i)}(t, \cdot)$ is regular outside the domain $D$, we conclude that $l_m^{(i)}(t; \cdot)$ is continuous on $H^1(\Omega)$ and satisfies

$$ |l_m^{(i)}(t; v)| \le C \|v\|_{H^1(\Omega)} $$

with a constant

$$ C = 2 \|\sigma\|_{L^\infty(\Omega)} \max_{t \in [0, T]} \|\nabla \tilde{u}_m^{(i)}(t, \cdot)\|_{L^2(\Omega \setminus \bar{D})} + C_\Gamma \sigma_D \max_{t \in [0, T]} \|\partial_n \tilde{u}_m^{(i)}(t, \cdot)\|_{L^2(\Gamma)} $$
which is independent from $t$, where $C_T > 0$ is the constant of the trace operator on $H^1(\Omega)$. The compatibility condition $l_m^{(i)}(t; 1) = 0$ can be proved as in [29] with the help of the solid angle formula. Then, Lax-Milgram's theorem guarantees existence and uniqueness of the solution at fixed time $t$. Now, notice that $\tilde{u}_m^{(i)}$ is continuous with respect to time on $\Omega \setminus D$ and $l_m^{(i)}$ depends continuously on $u_m^{(i)}$. This proves that $w_m^{(i)} \in C^0([0, T]; V)$ and we get

$$ \|w_m^{(i)}(t, \cdot)\|_V \le \frac{C}{\alpha} \quad \forall t \in [0, T] $$

where $\alpha > 0$ is the coercivity constant of $a(\cdot, \cdot)$.

<figure><img src="image_3.png" /></figure>

Now, the time-dependent EEG problem (2.14) with source term (2.13) admits a solution in the following sense:

**Theorem 2.2.** Let $\sigma \in L^\infty(\Omega)$ be such that $0 < \sigma_{\min} \le \sigma \le \sigma_{\max}$ a.e. on $\Omega$ with constants $\sigma_{\min}$ and $\sigma_{\max}$. Assume $(H_1) - (H_3)$ to be true. For $i \in \{1, \dots, d\}$ and $m \in \{1, \dots, M\}$, let $w_m^{(i)}$ be the solution of (2.21) where $\tilde{u}_m^{(i)}$ is defined by (2.18). Finally, denote by $q_m^{(i)}$ the $i$-th component of the moment $q_m \in L^2(0, T)^d$. Then,

$$ u = \sum_{m=1}^{M} \sum_{i=1}^{d} q_m^{(i)} (\tilde{u}_m^{(i)} + w_m^{(i)}) \quad (2.22) $$

is solution of the time-dependent EEG problem

$$ \begin{cases} \nabla \cdot (\sigma \nabla u) = F & \text{in } (0, T) \times \Omega, \\ \sigma \partial_n u = 0 & \text{on } (0, T) \times \Gamma, \end{cases} $$

with source term

$$ F(t, \cdot) = \sum_{m=1}^{M} q_m(t) \cdot \nabla \delta(\cdot - S_m(t)). $$

*Proof.* In order to prove that the potential $u$ under the form (2.22) is solution to problem (2.14), it is sufficient to notice that $\nabla \cdot (\sigma \nabla (\tilde{u}_m^{(i)} + w_m^{(i)}))$ is well defined at any time $t$ since

$$ \nabla \cdot (\sigma \nabla (\tilde{u}_m^{(i)} + w_m^{(i)})) = \nabla \cdot ((\sigma - \sigma_D)\nabla \tilde{u}_m^{(i)}) + \sigma_D \Delta \tilde{u}_m^{(i)} + \nabla \cdot (\sigma \nabla w_m^{(i)}) = F_m^{(i)}. $$

We further have $\sigma \partial_n (\tilde{u}_m^{(i)} + w_m^{(i)}) = 0$ on $[0, T] \times \Gamma$ by construction.

<figure><img src="image_4.png" /></figure>

### 3. DIFFUSIVE OPTICAL TOMOGRAPHY WITH TIME-DEPENDENT COEFFICIENTS

#### 3.1. Dimensional analysis of the P1-approximation of the radiative transfer equation

Light propagation is mainly governed by absorption and scattering phenomena. A popular model for photon transport in a context where the wavelength in the near infrared range is small (750–900 nm) compared to the characteristic length of the human head, is the radiative transfer equation (RTE).
Let $\mathbb{S}$ be the unit sphere in $\mathbb{R}^3$ and denote by $f$ a given phase function which describes the probability that a photon arriving from direction $\mathbf{s}$, is scattered in direction $\mathbf{s}'$. We assume that $f$ depends only from the angle between $\mathbf{s}$ and $\mathbf{s}'$ which is the case in an isotropic medium, and is such that

$$ \int_{\mathbb{S}} f(\mathbf{s} \cdot \mathbf{s}') \, d\mathbf{s}' = 1. $$

The unknown of the RTE is the specific intensity of light $L(t, \mathbf{x}, \mathbf{s})$ (expressed in $(\mathrm{W} \mathrm{m}^{-2} \mathrm{sr}^{-1})$) which depends on time $t$, position $\mathbf{x}$ and the direction $\mathbf{s}$. Let $Q$ be a light source. Then the RTE reads

$$ \frac{1}{c} \partial_t L(t, \mathbf{x}, \mathbf{s}) + \mathbf{s} \cdot \nabla L(t, \mathbf{x}, \mathbf{s}) + \mu_t L(t, \mathbf{x}, \mathbf{s}) = \mu_s \int_{\mathbb{S}} f(\mathbf{s} \cdot \mathbf{s}') L(t, \mathbf{x}, \mathbf{s}') \, d\mathbf{s}' + Q(t, \mathbf{x}, \mathbf{s}). \quad (3.1) $$

Here, $\mu_t = \mu_a + \mu_s$ is the total attenuation coefficient depending on the absorption coefficient $\mu_a$ and the scattering coefficient $\mu_s$ in the biological tissue, and $c$ is the speed of light in the medium.

Since the numerical resolution of the RTE results in a high computational cost, especially in the context of medical imaging applications, a common choice is to use the diffusion equation as an approximation. We will discuss here how to deal with time-dependent optical parameters. To this end, we define respectively, the diffuse photon density $\phi$ and the diffuse photon flux $\mathcal{L}$ by

$$ \phi(t, \mathbf{x}) = \int_{\mathbb{S}} L(t, \mathbf{x}, \mathbf{s}) \, d\mathbf{s}, \quad (3.2) $$

and

$$ \mathcal{L}(t, \mathbf{x}) = \int_{\mathbb{S}} L(t, \mathbf{x}, \mathbf{s}) \, \mathbf{s} \, d\mathbf{s}. \quad (3.3) $$

$\phi$ and $\mathcal{L}$ are respectively the moments of order 0 and 1 of the intensity $L$ which may be interpreted as the mathematical expectation and variance of $L$. We further assume that the source term $Q$ is isotropic, i.e. $Q$ does not depend on the direction $\mathbf{s}$ and thus $\int_{\mathbb{S}} Q(t, \mathbf{x}, \mathbf{s}) \, d\mathbf{s} = 0$. Expansion of the intensity $L$ in terms of spherical harmonics up to the first order then yields the following system

$$ \frac{1}{c} \partial_t \phi(t, \mathbf{x}) + \nabla \cdot \mathcal{L}(t, \mathbf{x}) + \mu_a(t, \mathbf{x}) \phi(t, \mathbf{x}) = Q \quad (3.4) $$

$$ \frac{1}{c} \partial_t \mathcal{L}(t, \mathbf{x}) + \frac{1}{3} \nabla \phi(t, \mathbf{x}) + (\mu_a(t, \mathbf{x}) + \mu_s'(t, \mathbf{x})) \mathcal{L}(t, \mathbf{x}) = 0, \quad (3.5) $$

where $\mu_s' = (1 - g)\mu_s$ is the reduced scattering coefficient, $g$ is the anisotropy factor and $Q(t, \mathbf{x}) = \int_{\mathbb{S}} Q(t, \mathbf{x}, \mathbf{s}) \, d\mathbf{s}$.

Since we are concerned with frequency modulated DOT systems, we consider a time-harmonic source term $Q(t, \mathbf{x}) = Q_0(\mathbf{x}) \exp(i\omega t)$ oscillating at a given pulsation $\omega > 0$. But the optical parameters $\mu_a$ and $\mu_s'$ depend both on the time and the space variable. We therefore cannot assume that the unknowns $\phi$ and $\mathcal{L}$ are separated in time and space as it is usually the case in time-harmonic problems. Let $\Phi(t, \mathbf{x})$ and $\mathbf{L}(t, \mathbf{x})$ be such that

$$ \phi(t, \mathbf{x}) = \Re(\Phi(t, \mathbf{x}) \exp(i\omega t)), $$

$$ \mathcal{L}(t, \mathbf{x}) = \Re(\mathbf{L}(t, \mathbf{x}) \exp(i\omega t)), $$

where $\Phi$ and $\mathbf{L}$ are complex-valued functions of time and space, and $i$ denotes the imaginary unit. System (3.4)-(3.5) then reads

$$ \frac{i\omega}{c} \Phi(t, \mathbf{x}) + \frac{1}{c} \partial_t \Phi(t, \mathbf{x}) + \nabla \cdot \mathbf{L}(t, \mathbf{x}) + \mu_a(t, \mathbf{x}) \Phi(t, \mathbf{x}) = Q_0(\mathbf{x}) \quad (3.6) $$

$$ \frac{i\omega}{c} \mathbf{L}(t, \mathbf{x}) + \frac{1}{c} \partial_t \mathbf{L}(t, \mathbf{x}) + \frac{1}{3} \nabla \Phi(t, \mathbf{x}) + (\mu_a(t, \mathbf{x}) + \mu_s'(t, \mathbf{x})) \mathbf{L}(t, \mathbf{x}) = 0. \quad (3.7) $$
<table><caption>TABLE 2. Baseline values for optical parameters of the neonatal head at 800 nm [33]. Constant anisotropy factor $g = 0.9$ and refractive index $n = 1.33$ for all tissues. Pulsation $\omega = 2\pi f$ with modulation frequency $f = 10^8$ (Hz). Homogenization of CSF/AT with $p = 15\%$ AT among CSF [36].</caption><thead><tr><th>Tissue</th><th>μ<sub>a</sub> (mm<sup>-1</sup>)</th><th>μ'<sub>s</sub> (mm<sup>-1</sup>)</th><th>μ<sub>a</sub> + μ'<sub>s</sub></th><th>ω/c (mm<sup>-1</sup>)</th></tr></thead><tbody><tr><td>scalp</td><td>0.018</td><td>1.9</td><td>1.918</td><td>0.0028</td></tr><tr><td>skull</td><td>0.016</td><td>1.6</td><td>1.616</td><td>0.0028</td></tr><tr><td>CSF/AT</td><td>0.0062</td><td>0.3122</td><td>0.3184</td><td>0.0028</td></tr><tr><td>grey matter</td><td>0.048</td><td>0.5</td><td>0.548</td><td>0.0028</td></tr><tr><td>white matter</td><td>0.037</td><td>1.0</td><td>1.037</td><td>0.0028</td></tr></tbody></table>

In the same way as for the time-dependent EEG problem, we now proceed to the scaling of system (3.6)–(3.7). Let $\tau$ and $\ell$ be the characteristic time scale and length of the problem. Let further $\Phi(t, \mathbf{x}) = \varphi\mathcal{P}'(t', \mathbf{x}')$ and $\mathbf{L}(t, \mathbf{x}) = \gamma\mathcal{L}'(t', \mathbf{x}')$ where $\varphi$ and $\gamma$ are reference quantities and $\mathcal{P}'$ and $\mathcal{L}'$ are dimensionless of order $\mathcal{O}(1)$ depending on $t' = t/\tau$ and $\mathbf{x}' = \mathbf{x}/\ell$. As before, prime notation corresponds to dimensionless operators or variables and system (3.6)–(3.7) becomes

$$ \frac{\omega}{c} \varphi \mathcal{P}' + \frac{1}{c\tau} \varphi \partial_{t'} \mathcal{P}' + \frac{\gamma}{\ell} \nabla' \cdot \mathcal{L}' + \mu_a \varphi \mathcal{P}' = Q_0, \quad (3.8) $$

$$ \frac{\omega}{c} \gamma \mathcal{L}' + \frac{1}{c\tau} \gamma \partial_{t'} \mathcal{L}' + \frac{1}{3\ell} \varphi \nabla' \mathcal{P}' + (\mu_a + \mu_s') \gamma \mathcal{L}' = 0. \quad (3.9) $$

Now, we take into account that frequency domain optical imaging devices employ radio-frequency modulated near-infrared light and operate at a typical frequency of $f = 10^8$ (Hz) [31]. According to [7, 32], oscillations in the concentrations of oxy- and deoxyhemoglobin occur at a frequency from 0.05 to 0.1 (Hz). Thus, two different time-scales appear simultaneously, a fast scale induced by the modulation frequency $f$ of the operating system and a slow scale $\tau$ resulting from neurovascular coupling. Since $\frac{1}{\tau} \ll \omega = 2\pi f$, dimensional analysis allows to justify that the second term in (3.8) (resp. (3.9)) is much smaller compared to the first term and we conclude that the time derivatives can be neglected.

We further investigate the orders of the first and last term in (3.9). Here below are listed typical baseline values for the optical parameters of the head tissues at a wavelength of 800 (nm). Except for the cerebrospinal fluid (CSF), these values are taken from [33]. The analysis of CSF has to be done with more caution. Indeed, it is commonly admitted that the diffusion approximation of RTE is much less precise in CSF than in the other tissues [34, 35]. This is mainly due to the small scale of the optical parameters which yields a large mean free path and consequently a large diffusion coefficient $\kappa$. However, it has been pointed out in [35, 36] that the subarachnoid space is not only filled with CSF but also contains thin filaments, the arachnoid trabeculae (AT), that connect the two inner meninges, arachnoid mater and pia mater. Taking into account these small scatterers by homogenization yields the values in Table 2 for a proportion of $p = 15\%$ AT among CSF (see [36] for details) and attests that $\omega/c \ll \mu_a + \mu_s'$ in all tissues. This implies that the first term in (3.9) can be neglected.

## 3.2. The time-dependent diffusion approximation

The system resulting from (3.8)–(3.9) and dimensional analysis reads

$$ \frac{\omega}{c} \Phi(t, \mathbf{x}) + \nabla \cdot \mathbf{L}(t, \mathbf{x}) + \mu_a(t, \mathbf{x}) \Phi(t, \mathbf{x}) = Q_0(\mathbf{x}), \quad (3.10) $$

$$ \frac{1}{3} \nabla \Phi(t, \mathbf{x}) + (\mu_a(t, \mathbf{x}) + \mu_s'(t, \mathbf{x})) \mathbf{L}(t, \mathbf{x}) = 0, \quad (3.11) $$
and implies that the photon flux $\mathbf{L}$ derives from the scalar potential $\Phi$. Substituting $\mathbf{L}$ in (3.10) then yields the diffusion equation

$$-\nabla \cdot (\kappa \nabla \Phi) + \left( \mu_a + \frac{\omega}{c} \right) \Phi = Q_0,$$

where the diffusion coefficient $\kappa(t, \mathbf{x}) = \frac{1}{3(\mu_a(t, \mathbf{x}) + \mu'_s(t, \mathbf{x}))}$ and the absorption coefficient $\mu_a(t, \mathbf{x})$ are functions of time and space variables. In order to take into account inner reflection at the boundary between the scalp and the exterior domain, a Robin-type boundary condition is prescribed on $\Gamma$ [16]. Hence, the boundary problem of the diffusion approximation with time-harmonic light source and time-dependent optical parameters reads

$$\left\{ \begin{array}{ll} -\nabla \cdot (\kappa \nabla \Phi) + \left( \mu_a + \frac{\omega}{c} \right) \Phi = Q_0 & \text{in } (0, T) \times \Omega, \\ \Phi + A\kappa \partial_n \Phi = 0 & \text{on } (0, T) \times \Gamma. \end{array} \right. \qquad (3.12)$$

Here, $A > 0$ is a scalar coefficient that depends on the refractive index (RI) mismatch between the biological tissue of RI $n$ and the surrounding medium of RI $n_a$. A precise formula for $A$ can be found in [31]. In comparison to the standard forward problem for DOT (e.g. [15, 37]), problem (3.12) is time-dependent and has been justified by dimensional analysis.

The variational formulation of problem (3.12) involves the complex vector spaces $L^2(\Omega; \mathbb{C})$ and $H^1(\Omega; \mathbb{C})$ of functions $v: \Omega \to \mathbb{C}$ defined on $\Omega$ with values in $\mathbb{C}$. In the sequel, $\bar{v}$ stands for the complex conjugate of $v$. Problem (3.12) is well posed in the following sense:

**Theorem 3.1.** Let $\mu_a$ and $\kappa$ belong to $C^0(0, T; L^\infty(\Omega))$ and assume that there are constants $\mu_{a,0} > 0$ and $\kappa_0 > 0$ such that

$$\mu_a \geq \mu_{a,0} \text{ and } \kappa \geq \kappa_0 \quad \forall t \in [0, T], \text{ a.e. on } \Omega. \qquad (3.13)$$

Let $Q_0 \in L^2(\Omega; \mathbb{C})$ and $A > 0$. Then, there is a unique solution $\Phi \in C^0(0, T; H^1(\Omega; \mathbb{C}))$ such that for any $t \in [0, T]$

$$a_t(\Phi(t, \cdot), v) = \ell(v) \quad \forall v \in H^1(\Omega; \mathbb{C}), \qquad (3.14)$$

where

$$a_t(u, v) = \int_{\Omega} \kappa(t, \cdot) \nabla u \cdot \nabla \bar{v} \, d\mathbf{x} + \int_{\Omega} \left( \mu_a(t, \cdot) + \frac{\omega}{c} \right) u \bar{v} \, d\mathbf{x} + \frac{1}{A} \int_{\Gamma} u \bar{v} \, ds \qquad (3.15)$$

is a sesqui-linear form defined for fixed $t$ on $H^1(\Omega; \mathbb{C}) \times H^1(\Omega; \mathbb{C})$ and

$$\ell(v) = \int_{\Omega} Q_0 \bar{v} \, d\mathbf{x} \qquad (3.16)$$

is an antilinear form on $H^1(\Omega; \mathbb{C})$.

*Proof.* Notice that the assumptions on the data $\mu_a, \kappa$ and $A$ guarantee that the sesqui-linear form $a_t(\cdot, \cdot)$ is continuous and coercive on $H^1(\Omega; \mathbb{C})$ for any $t \in [0, T]$ with constants independent from $t$. Since $Q_0 \in L^2(\Omega; \mathbb{C})$, the antilinear form $\ell$ is time-independent and continuous on $H^1(\Omega; \mathbb{C})$. At fixed time $t$, there is thus a unique solution $\Phi(t, \cdot) \in H^1(\Omega; \mathbb{C})$, solution to (3.14) due to Lax-Milgram's lemma for complex Hilbert spaces. Continuity of the solution $\Phi$ with respect to time then follows from the continuity of the parameters. Indeed,
let $(t_n)_{n \in \mathbb{N}}$ be a convergent sequence in $[0, T]$ such that $\lim_{n \to \infty} t_n = t^*$. Let $\Phi_n = \Phi(t_n, \cdot)$ and $\Phi^* = \Phi(t^*, \cdot)$ the corresponding solutions of (3.14) at $t_n$ and $t^*$, respectively. Then, for any $v \in H^1(\Omega; \mathbb{C})$, we have

$$ 
\begin{aligned}
a_{t_n}(\Phi_n - \Phi^*, v) &= a_{t_n}(\Phi_n, v) - a_{t^*}(\Phi^*, v) + a_{t^*}(\Phi^*, v) - a_{t_n}(\Phi^*, v) \\
&= \ell(v) - \ell(v) + a_{t^*}(\Phi^*, v) - a_{t_n}(\Phi^*, v) \\
&\longrightarrow 0
\end{aligned}
$$

as $n \to \infty$ since the coefficients $\kappa$ and $\mu_a$ are continuous with respect to $t$. Then, we infer from the coercivity of $a_t(\cdot, \cdot)$ that

$$ \|\Phi_n - \Phi^*\|_{1,\Omega} \le \frac{1}{\alpha} \Re (a_{t_n}(\Phi_n - \Phi^*, \Phi_n - \Phi^*)) \to 0 $$

as $n \to \infty$ since the coercivity constant $\alpha = \min(\kappa_0, \mu_{a,0})$ is time-independent. $\square$

<figure><img src="image_5.png" /></figure>

## 4. A MATHEMATICAL MODEL FOR NEUROVASCULAR COUPLING

Neurovascular coupling describes the hemodynamic response of the brain to the electrical activity of a group of neurons. Roughly speaking, the neuronal activity induces an increase of the cerebral blood flow in order to satisfy the increasing need of oxygen in the activated area. Changes in the concentrations of the involved chromophores oxyhemoglobin [HbO₂] and deoxyhemoglobin [Hb] imply changes in the optical parameters $\mu_a$ and $\mu'_s$. In this section, we summarize the main compartments of neurovascular coupling, i.e. a model for neuronal activity based on the work of [21] and the balloon model [18, 19] for the hemodynamic response. The balloon model has been developed initially to describe the BOLD signal of functional MRI [18], but the correlations between the BOLD signal and the concentration of deoxyhemoglobin allow to use it in diffuse optical imaging, too. The full model allows to get the input moment $q_m(t)$ in the source term (2.13) of the forward EEG problem (2.14) and the time-dependent absorption coefficient $\mu_a$ in the forward DOT problem (3.12) while explaining their link in neurovascular coupling. We start with the description of the neuronal model that yields the moment $q_m(t)$.

### 4.1. The neuronal model

Neurons are responsible for transmitting information between brain cells, through electrical impulses and chemical signals. The neuron is composed by three parts, the soma or cell body, the dendrites and the axon. When a stimulus occurs, the dendrites transmit the message to the soma which generates an action potential under appropriate conditions. The action potential is in turn sent through the axon and arrives at the presynaptic neuronal membrane which is separated from the postsynaptic neuron by the synaptic cleft. The action potential induces the release of neurotransmitters into the synaptic cleft which passes the information to the postsynaptic neuron and creates the postsynaptic potential.

Here, we present the neuronal model as a system of coupled ordinary differential equations. At the level of a single neuron, the membrane potentials $\nu_s$ and $\nu_d$ of the soma and the dendrite in presence of an input stimulus $I_{\text{stim}}$ are solutions of the following differential equations:

$$ C_m \frac{d\nu_s}{dt} = -I_{s,\text{tot}} + C_a(\nu_d - \nu_s) + I_{\text{stim}}, \quad (4.1) $$

$$ C_m \frac{d\nu_d}{dt} = -I_{d,\text{tot}} + C_a(\nu_s - \nu_d), \quad (4.2) $$

where $C_m$ is the specific membrane capacitance and $C_a$ is a constant which depends on the properties of the neuron.
Notice that we assume throughout the whole section that the quantities involved at the neuron level are spatially constant over the small area of the neuron membranes. This model could be generalized to spatially extended neurons by considering the cable equation or more general reaction-diffusion systems (see e.g. [38]), still resolved at the scale of a few neurons, resulting in much heavier computations.

## Voltage gated membrane currents and ionic concentrations

The total cross membrane currents $I_{*,tot}$ ($* \in \{s, d\}$) in the soma (resp. dendrite) compartment are the sum of voltage dependent sodium (Na) and potassium (K) currents $I_{s,Ion}$ (resp. $I_{d,Ion}$) and sodium, potassium and chlore leak currents as well as sodium-potassium exchange currents. Except for the leak currents, the ionic cross membrane currents are modelled using the Goldman-Hodgkin-Katz (GHK) formula

$$I_{*,Ion} = g_{*,Ion} \frac{F \nu_{*} ([Ion]_{i,*} - \exp(-\nu_{*}/\Phi) [Ion]_{e})}{\Phi \exp(-\nu_{*}/\Phi)}, \quad Ion \in \{Na, K\}, * \in \{s, d\}, \quad (4.3)$$

where $g_{*,Ion}$ is the channel specific conductance, $[Ion]_{i,*}$ and $[Ion]_{e}$ denote the (compartment specific) intra- and extracellular ionic concentrations, $F$ is the Faraday constant, and $\Phi = RT/F$ depends on the universal gas constant $R$ and the absolute temperature $T$.

The ionic exchange rates between the soma and the dendrite depend on the total ionic current and the flow of ions between the soma and the dendrite compartment which is assumed to be proportional to the difference of the ionic concentrations:

$$V_s \frac{d[Ion]_{i,s}}{dt} = -C_s I_{s,Ion} + C_{Ion} ([Ion]_{i,d} - [Ion]_{i,s}) \quad (4.4)$$

$$V_d \frac{d[Ion]_{i,d}}{dt} = -C_d I_{d,Ion} + C_{Ion} ([Ion]_{i,s} - [Ion]_{i,d}) \quad (4.5)$$

$$\frac{V_e}{V_s + V_d} \frac{d[Ion]_e}{dt} = \frac{C_s}{V_s} I_{s,Ion} + \frac{C_d}{V_d} I_{d,Ion} \quad (4.6)$$

where $V_s$, $V_d$ and $V_e$ are respectively the volume of the soma, the dendrite and the extracellular space, and $C_s$, $C_d$ and $C_{Ion}$ are constants depending on the neuron and the ion under consideration.

The ionic currents given by the GHK-formula (4.3) travel across voltage gated ionic channels. Each channel has a conductance $g_{*,Ion} = m_{*,Ion}^p h_{*,Ion}^q \bar{g}_{Ion}$, product of an ion-specific maximum conductance $\bar{g}_{Ion}$ and a factor $m_{*,Ion}^p h_{*,Ion}^q \in [0, 1]$ which describes the fraction of active and inactive open channels. The variables $m_{*,Ion}$ and $h_{*,Ion}$ are solutions of a couple of differential equations

$$\frac{dm_{*,Ion}}{dt} = \alpha_{m,*,Ion}(\nu_*)(1 - m_{*,Ion}) - \beta_{m,*,Ion}(\nu_*)m_{*,Ion}, \quad (4.7)$$

$$\frac{dh_{*,Ion}}{dt} = \alpha_{h,*,Ion}(\nu_*)(1 - h_{*,Ion}) - \beta_{h,*,Ion}(\nu_*)h_{*,Ion}, \quad (4.8)$$

with opening and closing rates $\alpha$ and $\beta$ depending on the membrane potential $\nu_*$.

## Leak currents

The leak currents are given by Ohm's law,

$$I_{*,leak,Ion} = \bar{g}_{Ion}(\nu_{*} - E_{Ion}) \quad (4.9)$$

where $E_{Ion}$ denotes the ion-specific potential at rest.
## Sodium-potassium exchange currents

The coupling between the neuronal and the vascular compartment is realized by the sodium-potassium exchange pump which aims to balance the ionic concentrations between the intra- and extracellular spaces. These pumps are present as well in the soma as in the dendrite compartment. The pump moves out three intracellular sodium ions and two extracellular potassium ions:

$$I_{*,Na,pump} = 3I_{*,pump}, \qquad (4.10)$$

$$I_{*,K,pump} = -2I_{*,pump}, \qquad (4.11)$$

where the total current $I_{*,pump}$ due to the Na/K-exchange pump is given as the product of a function $\gamma_{1,*}$ depending on the concentrations of intracellular sodium and extracellular potassium and a function $\gamma_{2,*}$ depending on the tissue oxygen concentration $[O_2]$,

$$I_{*,pump} = I_{\max}\gamma_{1,*}([K]_e, [Na]_{*,i})\gamma_{2,*}([O_2]). \qquad (4.12)$$

Expressions for $\gamma_{1,*}$ and $\gamma_{2,*}$ can be found in [21].

In turn, the variation of oxygen concentration $[O_2]$ depends on the oxygen supply and consumption, and satisfies the differential equation

$$\frac{d[O_2]}{dt} = J \frac{[O_2]_b - [O_2]}{[O_2]_b - [O_2]_0} - M. \qquad (4.13)$$

Here, $J$ is the rate of change in oxygen concentration, $[O_2]_b$ and $[O_2]_0$ are, respectively, the steady state and baseline oxygen concentration, and $M = M([O_2], [K]_e, [Na]_{i,s}, [Na]_{i,d})$ denotes the cerebral metabolic rate of oxygen (CMRO$_2$) which accounts for variation in time of the oxygen consumption after stimulation due to the exchange pump and other processes [21].

## 4.2. The vascular compartment

One of the most popular models to describe the dynamics of the biomechanical and metabolic quantities involved in brain activation is the balloon model [18, 19]. It consists of two coupled differential equations for the concentration of deoxygenated hemoglobin [Hb] and the cerebral blood volume CBV of the venous compartment:

$$\frac{d[\text{Hb}]}{dt} = \frac{1}{\tau_0} \left( \frac{M}{M_0} - \frac{[\text{Hb}]}{\text{CBV}} f_{\text{out}} \right), \qquad (4.14)$$

$$\frac{d\text{CBV}}{dt} = \frac{1}{\tau_0} \left( \frac{\text{CBF}}{\text{CBF}_0} - f_{\text{out}} \right). \qquad (4.15)$$

Here, $\tau_0$ is the mean transit time at rest, i.e. the time that a given blood volume stays in the capillary circulation, and $M_0$ and $CBF_0$ are equilibrium values for the metabolic rate of oxygen $M$ and the cerebral blood flow $CBF$, respectively. The cerebral blood flow in the arterioles of the vascular compartment is of Poiseuille type,

$$\text{CBF} = \text{CBF}_0 \frac{R^4}{R_0^4} \qquad (4.16)$$

where $R$ is the radius of the arterioles and $R_0$ is the corresponding steady state value. The outgoing flow $f_{\text{out}}$ is related to the blood volume CBV by

$$f_{\text{out}} = \text{CBV}^{\frac{1}{d}} + \tau \frac{d\text{CBV}}{dt} \qquad (4.17)$$

with empirical constants $d$ and $\tau$.
### 4.3. The complete differential system

In summary, we get a system of ordinary differential equations

$$ \frac{dU}{dt} = \mathcal{F}(U) \qquad (4.18) $$

which has to be completed by initial values $U(0) = U_0$. The vector $U$ of unknowns can be split into two parts, on the one hand the neuronal part $U_{\text{neuron}}$ which contains the membrane potentials $\nu_s$ and $\nu_d$, the intra- and extracellular ionic concentrations $[\text{Ion}]_{i,*}$ and $[\text{Ion}]_e$ and finally the gating variables $m$ and $h$, on the other hand the hemodynamic part $U_{\text{hd}}$ containing the concentrations of oxygen and deoxygenated hemoglobin $[\text{O}_2]$ and $[\text{Hb}]$ as well as the cerebral blood volume CBV. The right hand side of (4.18) may also be split into a vector $\mathcal{F}_{\text{neuron}}$ and $\mathcal{F}_{\text{hd}}$ but due to neurovascular coupling, $\mathcal{F}_{\text{neuron}}$ acts on $U_{\text{neuron}}$ as well as on $[\text{O}_2]$ through the presence of the exchange pump current $I_{*,\text{pump}}$, whereas $\mathcal{F}_{\text{hd}}$ acts on $U_{\text{neuron}}$ as well as on $U_{\text{hd}}$ since the metabolic rate of oxygen $M$ depends on the ionic concentrations $[\text{K}]_e$ and $[\text{Na}]_{i,*}$. Thus, system (4.18) reads

$$ \frac{d}{dt} \begin{pmatrix} U_{\text{neuron}} \\ U_{\text{hd}} \end{pmatrix} = \begin{pmatrix} \mathcal{F}_{\text{neuron}}(U_{\text{neuron}}, [\text{O}_2]) \\ \mathcal{F}_{\text{hd}}(U_{\text{neuron}}, U_{\text{hd}}) \end{pmatrix}. \qquad (4.19) $$

Notice that the full system of differential equations described in [20, 21] contains other compartments than the two mentioned here, as for example astrocyte.

### 4.4. The time-dependent source term of the forward EEG-problem

#### The post-synaptic current

In response to an action potential $\nu_s$ (i.e. a presynaptic spike), neurotransmitters are released at the presynaptic terminal. These neurotransmitters then move into the synaptic cleft and bind with receptors in the postsynaptic neuron, opening ion channels. This movement of ions across the neuronal membrane generates a postsynaptic current $I_{\text{post}}$ which is given by

$$ I_{\text{post}}(t) = g(t)(\nu_{\text{post}}(t) - \nu_{\text{rev}}). \qquad (4.20) $$

Here, the conductance $g(t)$ depends on the presynaptic neuron, $\nu_{\text{post}}$ is the potential across the postsynaptic membrane, and $\nu_{\text{rev}}$ is a (constant) reversal potential the value of which depends on the neurotransmitter.

There are different types of neurotransmitters and the model will depend on which neurotransmitter we consider. Here, we focus on a single channel model which is valid e.g. for AMPA or GABA$_A$ which are, respectively, excitatory or inhibitory neurotransmitters. As before, we express the conductance as a maximum conductance $\bar{g}$ and a factor between 0 and 1 modelling the proportion of open channels. Thus, we have $g(t) = \bar{g}s(t)$ in (4.20) where $s$ is obtained by solving the differential equation [41]

$$ \frac{ds}{dt} = K_1[T](1-s) - K_2s, \qquad (4.21) $$

with $K_1, K_2$ two constants that depend on the neurotransmitter and $[T]$ the concentration of neurotransmitters in the synaptic cleft. We use a relationship between $[T]$ and $\nu_s$ from [41] according to which $[T]$ is given by

$$ [T](\nu_s) = \frac{T_{\text{max}}}{1 + \exp(-(\nu_s - V_T)/K_p)}, \qquad (4.22) $$
where $T_{max}$ is the maximum concentration of transmitters in the synaptic cleft, $V_T$ the value at which the concentration is halfed, $K_p$ models the steepness and $\nu_s$ is the presynaptic action potential, given by the resolution of (4.18).

Finally, the corresponding postsynaptic current is then given by the right hand side of (4.20) where the postsynaptic potential is taken as the membrane potential $\nu_d$ in the dendrites.

### The time-dependent current dipolar moment

The above analysis describes potentials and currents for a single neuron. EEG signals at the scalp, however, result from postsynaptic currents generated by a large amount of synchronized pyramidal neurons with similar orientations which sum up according to Kirchhoff's law.

The source term of the time-dependent forward EEG problem (2.14) is the divergence of the source current density $\mathbf{j}^p$ which we assume here to result from one electrical current dipole, i.e. $M = 1$. If the active area of the brain is time-independent and localized at the point $S$, this implies, according to Section 2, that $\mathbf{j}^p$ is separated in the variables $t$ and $\mathbf{x}$ and may be written as a dipolar moment $\mathbf{q}$ depending on the time variable $t$, multiplied by a Dirac delta $\delta_S$, depending only on $\mathbf{x}$. The moment $\mathbf{q}$ depends on the postsynaptic current $I_{post}$, a length scale $\Delta s$ and the number $N_{syn}$ of synchronized neurons,

$$ \mathbf{j}^p = \mathbf{q}(t)\delta_S = N_{syn}I_{post}(t)\Delta s\mathbf{e}(t)\delta_S, \quad (4.23) $$

where $\mathbf{e}(t) \in \mathbb{R}^3$ is a unit vector describing the orientation of the neurons at time $t$. The solution $u$ of (2.14) then models the spatial distribution in the head of the electric potential generated by the neuronal current $I_{post}$ at $S$ over a time interval $(0, T)$.

The order of magnitude of $\mathbf{j}^p$ is of 10 (nA m$^{-2}$) for the following data

$$ N_{syn} = 10^6, \quad I_{post} \approx 10(\text{pA}), \quad \Delta s = 1(\text{mm}). $$

The corresponding order of magnitude of the simulated potential $u$ at the electrodes is about 10 (µV) which is consistent with experimental data (cf. e.g. [42]).

## 4.5. Time-dependent optical parameters

Contrary to [21], we are not interested here in the BOLD signal, but in the concentrations of deoxy- and oxyhemoglobin [Hb]$_N$ and [HbO$_2$]$_N$ which we assume normalized by their (initial) values at rest. [Hb] is part of the vector of unknowns of system (4.18). [HbO$_2$]$_N$ can be obtained by the formula

$$ [\text{HbO}_2]_N = [\text{HbT}]_N - [\text{Hb}]_N + 1 \quad (4.24) $$

where the total normalized concentration of hemoglobin [HbT]$_N$ is given by

$$ [\text{HbT}]_N = \frac{[\text{Hb}]_N \text{CBF}_N}{M_N}, \quad (4.25) $$

with normalized cerebral blood flow CBF$_N$ and normalized change of oxygen consumption $M_N$.

The main chromophores contributing to absorption in the biological tissues of the brain are deoxy- and oxyhemoglobin as well as water. The absolute absorption coefficient at rest can therefore be related to the concentrations of the chromophores by

$$ \mu_{a,0} = \mu_{a,\text{water}} P_{\text{water}} + \alpha_{\text{Hb}} [\text{Hb}]_0 + \alpha_{\text{HbO}_2} [\text{HbO}_2]_0, \quad (4.26) $$
<figure><img src="image_6.png" /><figcaption>FIGURE 2. The full pipeline for simulation of coregistered EEG- and DOT-signals: a simplified representation of the different compartments of the ODE system (Neuronal compartment (4.1)–(4.9), Vascular compartment (4.14)–(4.17), Coupling (4.10)–(4.13)), Source term (divergence of term (4.23)), $\Delta\mu_a$ (4.27), EEG signal (2.18)–(2.21)–(2.22), DOT signal (3.14).</figcaption></figure>

<table><caption>TABLE 3. Characteristics of the realistic head mesh.</caption><thead><tr><th>Nodes</th><th>Tetrahedra</th><th>h<sub>min</sub> (mm)</th><th>h<sub>max</sub> (mm)</th><th>diameter (mm)</th></tr></thead><tbody><tr><td>108 669</td><td>590 878</td><td>0.34</td><td>14</td><td>120</td></tr></tbody></table>

where $P_{water}$ is the proportion of water in the cortex, and $\alpha_{Hb}$ and $\alpha_{HbO_2}$ are the specific extinction coefficients of Hb and $HbO_2$ [43]. Changes in deoxy- and oxyhemoglobin concentrations during cerebral activation lead to the modification of $\mu_{a,0}$ which is expressed in terms of the time-dependent normalized concentrations $[Hb]_N$ and $[HbO_2]_N$:

$$ \Delta\mu_a(t) = \alpha_{Hb}[Hb]_0([Hb]_N(t) - 1) + \alpha_{HbO_2}[HbO_2]_0([HbO_2]_N(t) - 1). \quad (4.27) $$

The effect of cerebral activity on the scattering coefficient $\mu_s$ is less documented and will be neglected here.

## 5. NUMERICAL SIMULATION OF COREGISTERED EEG/DOT SIGNALS

In this section, we explain the numerical pipeline from an input stimulation current up to simulated coregistered EEG/DOT signals at the electropotodes according to Figure 2. Electropotodes are a medical device which allows to monitor simultaneously the neuronal and vascular activities of the brain [6].

We provide numerical simulations on a realistic head model of a healthy fullterm newborn obtained from coregistration of MR and CT images of the Amiens' hospital database (courtesy GRAMFC, INSERM U1105, Amiens, France (H. Azizollahi [27])). The model distinguishes between four tissues: brain, CSF, skull, and scalp. Table 3 summarizes the mesh parameters. Tables 1 and 2 give respectively the conductivity $\sigma$ and the optical parameters $\mu_{a,0}$ and $\mu'_s$ of the tissues.

Synthetic measurements are observed at 8 pointwise electropotodes $\{e_l\}_{l=1}^8$ situated in different regions of the scalp on the boundary $\Gamma$. In Figure 3 (left), one distinguishes the electropotodes $e_3$ and $e_5$ situated respectively in the frontal and central area as well as $e_2$ and $e_8$ which are placed above the left hemisphere. We numerically simulate the measured EEG/DOT signals at electropotodes during the observation time interval $[0, T]$, with $T = 150$ (s). The ODE system yields information about the time dependence of the moment $q(t)$ and the changes in the optical parameter $\Delta\mu_a(t)$. The position $S$ of the electrical brain activity as well as the spatial distribution of $\mu_a(t, x)$ with respect to the value at rest $\mu_{a,0}$ have to be fixed ad hoc. We consider a single source point $S$ below the electropotde $e_2$ and a moment $q(t)$ directed normally to the brain/LCS interface (Fig. 3 middle and right). We assume that the changes in the absorption parameter occur in a ball $B_r(S)$ of radius $r = 12$ (mm) which is centered at the point $S$ of electrical activity (see Fig. 4). The spatial distribution of $\Delta\mu_a$
<figure><img src="image_7.png" /><figcaption>FIGURE 3. Position of the electroptodes (left). Frontal and sagittal view of position (black point) and orientation of the source.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 4. Localization in the brain of the perturbation in the absorption coefficient. Support of $\Delta\mu_a$ in red. $\Delta\mu_a$ is constant in the support.</figcaption></figure>

inside the perturbation is assumed to be constant such that $\mu_a(t, \mathbf{x})$ is given by the formula

$$ \mu_a(t, \mathbf{x}) = \mu_{a,0}(t, \mathbf{x}) + \Delta\mu_a(t)\chi_{B_r(S)}(\mathbf{x}), \quad (5.1) $$

where $\chi$ is the characteristic function of the ball $B_r(S)$. The volume of the ball is comparable to the one in [33] and represents 1.6% of that of the brain. It models a small and localized perturbation of the absorption coefficient caused by changes in deoxy- and oxyhemoglobin concentrations in the area of the electrical activity. *Numerical pipeline*. We provide detailed information on the implementation of each part of the method.

**Step 1:** ODE system. The method consists first in solving the complete system (4.18) of differential equations which, recall, models the neurovascular coupling at the neuron level (Neuronal compartment + Vascular compartment + Coupling). To this end, we apply the ode15s Matlab solver, which is an implicit variable-step variable-order method (see [44] for details on the involved numerical differentiation formulas) intended for stiff problems. Indeed, we need to simulate phenomena that occur at different time scales. The neuronal activity takes a few milliseconds and can show strong variations. The vascular part is characterized by a smooth behaviour along a period of several seconds. An order of 10000 time steps ($t_k$)$_{1 \le k \le K}$ is required to obtain accurate results.
Then, the EEG and DOT problems are solved simultaneously on the head mesh (see Tab. 3) during the time interval $[0, T]$ with input data obtained from the solution of system (4.18).

Step 2: EEG signals. We compute the postsynaptic current $I_{\text{post}}(t_k)$ from (4.20) (see Sect. 4.4) and deduce the time-dependent source term $F = \nabla \cdot \mathbf{j}^p$ of the forward EEG model (2.14) from formula (4.23). We solve problem (2.14) by the subtraction approach and use the finite element method with Lagrange-$\mathbb{P}_1$ elements implemented in FreeFem++ [45] to compute the regular potentials $w_m^{(i)}$, $i \in \{1, \dots, d\}$, solution to the variational formulation (2.21). Notice, that the functions $w_m^{(i)}$ are time-independent in the present context of a static source position $S$. They require the resolution of $d$ linear systems, $d = 3$ in the present three-dimensional configuration, involving the same matrix but different right hand sides. The linear systems are solved by a direct solver for sparse matrices. This allows to compute only once the underlying factorization of the system's matrix. The potential $u$ is then computed by linear combination from (2.22) and EEG measurements are given by

$$U_{k,l} = u(t_k, \mathbf{e}_l), \quad k = 1, \dots, K, \quad l = 1, \dots, 8, \qquad (5.2)$$

where $K$ is the number of time steps. Due to the Neumann boundary condition, the solution $u$ is determined up to an additive constant only. In Section 2, we considered the subspace of functions in $H^1(\Omega)$ with vanishing mean value to get uniqueness of the solution. In practice, the absolute potential cannot be measured by EEG and the measured voltage is the difference between the potential at electropotode $\mathbf{e}_l$ and the potential at a reference electropotode. Numerically, we therefore impose a Dirichlet boundary condition on a single node of the mesh to fix the constant. Since the singularity of the potential $u$ is taken into account exactly in the subtraction approach, the discretization error is given by the error on the regular potential. We refer to [29] for a detailed error analysis.

Step 3: DOT signals. Concerning the vascular part, the resolution of system (4.18) allows the computation of the variation $\Delta\mu_a$ at time $t_k$ from formula (4.27) (see Sect. 4.5) and the perturbed absorption coefficient is given by

$$\mu_a(t_k, \mathbf{x}) = \mu_{a,0}(t_k, \mathbf{x}) + \Delta\mu_a(t_k)\chi_{B_r(S)}(\mathbf{x}),$$

according to (5.1). The light source $Q_0$ is modeled as a Gaussian function

$$Q_0(\mathbf{x}) = A_s \exp \left( -\frac{\|\mathbf{x} - \mathbf{x}_s\|}{2\sigma_s^2} \right)$$

where $\sigma_s = 5 \times 10^{-3}$ fixes the beam diameter and the amplitude is set to $A_s = (2\pi\sigma_s)^{-1}$. $\mathbf{x}_s$ is the center of the source which should be located at a distance $\delta \sim 1/\mu'_s$ away from the surface $\Gamma$ in order to justify an isotropic source term. We use Lagrange-$\mathbb{P}_1$ finite elements implemented in FreeFem++ to compute an approximation of the solution to the variational formulation (3.14). Notice that (3.14) requires the assembly of the matrix associated with the sesqui-linear form $a_t(\cdot, \cdot)$ as well as the resolution of the underlying linear system at any discrete time step $t_k$. We employ a coarsened time grid in order to reduce the computational time. This is numerically justified since variation in time for DOT arises at a larger time scale than in EEG. We normalize the photon density $\phi$ with respect to the reference configuration with no perturbation, i.e $\mu_a(t_k, \mathbf{x}) = \mu_{a,0}$. DOT measurements are thus expressed by

$$\phi_{k,l} = \frac{|\phi(t_k, \mathbf{e}_l)|}{|\phi_0(t_k, \mathbf{e}_l)|}, \quad k = 1, \dots, K', \quad l = 1, \dots, 8,$$

where $\phi_0$ is the solution of the DOT problem (3.14) with piecewise constant absorption coefficient $\mu_{a,0}$ and $K' \ll K$ is the number of reduced time steps. Baseline values for $\mu_{a,0}$ and $\mu'_s$ are those from Table 2.
<figure><img src="image_9.png" /><figcaption>FIGURE 5. (A) Simulated postsynaptic current $I_{\text{post}}$ for a rectangular pulse stimulation of amplitude $I_{\text{stim,max}} = 22$ (mA) between 100 and 101 (s). (B) Zoom of $I_{\text{post}}$ to the interval [100, 101]. (C) Normalized concentrations in deoxy- and oxyhemoglobin for the same stimulation.</figcaption></figure>

*Simulations at the neuron level.* We first analyse the link between the input stimulus $I_{\text{stim}}$ and the output of the postsynaptic current $I_{\text{post}}$ on the one hand and variations in the normalized concentrations [Hb]$_N$ and [HbO$_2$]$_N$ of deoxy- and oxyhemoglobin on the other. The model parameters are those of the neurotransmitter GABA$_A$. According to [20, 21], [Hb]$_N$ has been obtained from [Hb] by normalization with respect to the mean of maximum and minimum values of [Hb] in the time interval $[T_1 - 20, T_1]$ where neural stimulation starts at time $T_1$. In Figure 5, we present the coupling between $I_{\text{post}}$ and the variations of [Hb]$_N$ and [HbO$_2$]$_N$ for a rectangular pulse stimulation $I_{\text{stim}}$ with amplitude $I_{\text{stim,max}} = 22$ (mA) which occurs during 1s in the time interval [100, 101]. We may notice that the most significant electrical activity is limited to the stimulation interval whereas an increase in [HbO$_2$]$_N$ and a symmetrical decrease in [Hb]$_N$ arise a few seconds after the neural activity.

We next study the influence of the amplitude of $I_{\text{stim}}$ on the postsynaptic current and the hemodynamic concentrations. In Figure 6, we present $I_{\text{post}}$ and the variations of [Hb]$_N$ and [HbO$_2$]$_N$ for a rectangular pulse stimulation $I_{\text{stim}}$ with different amplitudes, $I_{\text{stim,max}} = 22$ (mA) (left column) and $I_{\text{stim,max}} = 19$ (mA) (middle). We notice that the amplitude of the stimulation pulse influences the spike density of the postsynaptic current as it was mentioned in [21, 46]. This observation can be related to results in [41, 47] where the authors state that the generation of an action potential shows an all-or-none behaviour in the sense that below a given threshold of the applied current amplitude, no action potential will occur. If the applied current is strong enough during a sufficiently long time interval, a periodic action potential (a "burst of spikes") is generated. According to [41], these spikes arise also in both the postsynaptic potential and postsynaptic current at the neuron level. Results in the right column of Figure 6 correspond to an input stimulus $I_{\text{stim}}$ containing two rectangular pulses (the first one of amplitude $I_{\text{stim,max}} = 22$ (mA) and the second $I_{\text{stim,max}} = 19$ (mA)). The postsynaptic current has two periods of spikes, associated to the two pulses, with a higher density for the first one. One may notice that the global behaviour of the hemodynamic response (lower line of Fig. 6) is similar for the three configurations even if the concentration of oxyhemoglobin shows two nearby peaks corresponding to the two periods of spikes in the postsynaptic current.
<figure><img src="image_10.png" /><figcaption>FIGURE 6. Simulated postsynaptic currents $I_{\text{post}}$ (middle line, zoom on interval [100, 101]) and normalized concentrations of deoxy- and oxyhemoglobin (lower line) for different stimulation currents (upper line): rectangular pulse with amplitude $I_{\text{stim,max}} = 22$ (mA) (left) and $I_{\text{stim,max}} = 19$ (mA) (middle), two rectangular pulses of amplitudes 22 (mA) and 19 (mA) (right).</figcaption></figure>

*Simulations of EEG-DOT signals.* We have numerically exhibited the neurovascular coupling at the neuron level. Now, we are going to simulate the corresponding EEG/DOT signals at the electropodes. We consider the rectangular pulse stimulation $I_{\text{stim}}$ of amplitude $I_{\text{stim,max}} = 22$ (mA). In Figure 7, we report two views of the electric potential $u(t_k, \cdot)$, solution to the EEG problem at time $t_k = 108$ s. Figures 8 and 9 show the hemodynamic counterpart, i.e. the modulus of the photon density $\phi(t_k, \cdot)$, solution to the DOT problem. In Figure 10, we report the corresponding simulated EEG/DOT signals. One may notice that the magnitude of the simulated DOT signals depends linearly on the source term which should be calibrated precisely in case of comparison with experimental measurements. This point has not been investigated in the present study. It is also known that the amplitude of the DOT signal decreases rapidly with the distance to the source point which is confirmed by our simulation. A comparison between simulated signals and experimental measurements can be found in [48] for a cylindrical phantom with optical coefficients which are about three times smaller than those of scalp or skull. The reported normalized values vary between 1e-10 and 1.

According to formula (2.22), the electric potential $u$ depends linearly on $I_{\text{post}}$. Thus, EEG signals have the same temporal behavior than $I_{\text{post}}$ (see Fig. 5). Spikes are visible between $t = 100$ (s) and $t = 101$ (s) whereas the amplitude differs with respect to the position of the electrode. Electropode $e_2$, which is located near the brain source $S$, records the most important changes in both the electric potential $u$ and the photon density $\phi$ (between 0.15% and 15% according to Fig. 10). We observe that the simulated EEG and DOT signals are correlated in time. The numerical results of Figures 6 and 10 illustrate the neurovascular coupling in the sense that variation of the chromophores [Hb] and [HbO$_2$] starts at the onset of the EEG burst and takes several seconds before coming back to baseline values. This is consistent with the clinical observations in [32] where simultaneous EEG and NIRS recordings on premature infants during quiet sleep are discussed. One may notice
<figure><img src="image_11.png" /><figcaption>FIGURE 7. Frontal and sagittal views of the electric potential $u$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 8. Sagittal views of the source term (below e8) and the modulus $|\Phi|$ of the photon density (at e2).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 9. Frontal and sagittal views of the modulus $|\Phi|$ of the photon density.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 10. Numerical simulation of coregistered EEG and DOT signals at electrotodes $e_2$ (left), $e_3$ (middle) and $e_5$ (right) for a rectangular pulse stimulation of $I_{stim,max} = 22$ (mA). Upper line: electric potential $u$ (zoom on interval [100, 101]). Lower line: modulus of photon density $|\Phi|$ (in $\log_{10}$-scale).</figcaption></figure>

that [32] mentioned a beginning of the [Hb] and [HbO₂] variations 3 (s) before the EEG burst onset that we do not observe in the simulations. The recorded EEG traces in [32] correspond to events that occur without exterior stimulation and differ from our simulated ones. They show however a precise onset of electrical brain activity which makes the comparison meaningful. Another clinical study [7] related to spike-and-wave discharge during absence epilepsy shows EEG channels with a sharp onset and sharp end during the seizure as in our simulation. The corresponding [Hb] and [HbO₂] variations are however inverted in the sense that hemoglobin oxygenation increases several seconds before the onset. This observation has been made for epilepsy seizures wheras our simulation is not related to a particular event. We also cite a clinical study [49] which establishes a classification of hemodynamic responses related to EEG bursts for preterm neonates. It is noticed that variations in [Hb] and [HbO₂] can differ in shape from one patient to another: [Hb] may be increasing at the onset, [HbO₂] decreasing, or vice versa, or [Hb] and [HbO₂] are simultaneously increasing or decreasing.

## 6. CONCLUSION

In this paper, we have developed and analysed a mathematical and numerical pipeline for simulating coregistered EEG and DOT measurements. It combines several compartments: a non-linear system of differential equations (existing in the literature) which models the neurovascular coupling at the neuron level, and PDE-based derived from the standard stationary models for EEG and DOT that allow to compute time-dependent signals at the electrotodes on the scalp. Evolution in time has been introduced into the (stationary) EEG and DOT models by considering time-dependent source terms in EEG and time-dependent optical parameters for DOT. Dimensional analysis proved that time-derivatives can be neglected in the time-dependent PDEs. In particular, this shows that time-dependent sources (resp. parameters) do not modify the nature of the models, i.e. the stationary equations can be written over a time interval. We further assumed that the source term in the EEG model as well as the optical parameters of the DOT model are separated in time and space variables. This allows decoupling of the ODE system and the PDE models. For EEG, it corresponds to the common choice of an electrical current dipole. Existence and uniqueness results in appropriate vector spaces attest the well-posedness of both PDE models. We have detailed the pipeline from the input stimulation of cerebral activity up to the coregistered EEG/DOT signals. Numerical results on a three-dimensional realistic head model illustrate the electric potential and hemodynamic changes occurring from neurovascular coupling. Up to the best of our knowledge, this is the first contribution for modelling coupled EEG and DOT signals by solving both ODEs and PDEs. The generalization to other ionic models that provide both the postsynaptic current and the concentrations of deoxy- and oxyhemoglobin is straightforward. Neural models based e.g. on reaction-diffusion PDEs or
multi-compartmental approaches as well as spatio-temporal models for the concentration of the hemodynamic parameters, and thus of the absorption coefficient, are an interesting perspective for more involved coupled problems.

Coregistered EEG/DOT signals give valuable information to clinicians for a deepful understanding of brain disorders. The derivation of a EEG/DOT forward model is therefore an essential preliminary step to the resolution of the corresponding inverse problem which consists in reconstructing simultaneously the location of electric sources and the optical parameters in the brain.

DATA AVAILABILITY STATEMENT

No new data/codes were created or analyzed in this study.

REFERENCES

[1] J.M. Rennie, L.S. de Vries, M. Blennow, A. Foran, D.K. Shah, V. Livingstone, A.C. van Huffelen, S.R. Mathieson, E. Pavlidis, L.C. Weeke, M.C. Toet, M. Finder, R.M. Pinnamaneni, D.M. Murray, A.C. Ryan, W.P. Marnane and G.B. Boylan, Characterisation of neonatal seizures and their treatment using continuous EEG monitoring: a multicentre experience. *Arch. Dis. Child Fetal Neonatal Ed.* **104** (2019) F493–F501.

[2] T. Koskela, G.S. Kendall, S. Memon, M. Sokolska, T. Mabuza, A. Huertas-Ceballos, S. Mitra, N.J. Robertson, J. Meek and K. Whitehead, Prognostic value of neonatal EEG following therapeutic hypothermia in survivors of hypoxic-ischemic encephalopathy. *Clin. Neurophysiol.* **132** (2021) 2091–2100.

[3] C. Lee, R. Cooper and T. Austin, Diffuse optical tomography to investigate the newborn brain. *Pediatr. Res.* **82** (2017) 376–386.

[4] M. Nourhashemi, M. Mahmoudzadeh, S. Goudjil, G. Kongolo and F. Wallois, Neurovascular coupling in the developing neonatal brain at rest. *Hum. Brain Mapp.* **41** (2020) 503–519.

[5] A. Gallagher, F. Wallois and H. Obrig, Functional near-infrared spectroscopy in pediatric clinical research: different pathophysiologies and promising clinical applications. *Neurophotonics* **10** (2023) 023517.

[6] F. Wallois, M. Mahmoudzadeh, A. Patil and R. Grebe, Usefulness of simultaneous EEG-NIRS recording in language studies. *Brain Lang.* **121** (2012) 110–123.

[7] N. Roche-Labarbe, B. Zaaimi, P. Berquin, A. Nehlig, R. Grebe and F. Wallois, NIRS-measured oxy- and deoxyhemoglobin changes associated with EEG spike-and-wave discharges in children. *Epilepsia* **49** (2008) 1871–1880.

[8] M. Chalia, C.W. Lee, L.A. Dempsey, A.D. Edwards, H. Singh, A.W. Michell, N.L. Everdell, R.W. Hill, J.C. Hebden, T. Austin and R.J. Cooper, Hemodynamic response to burst-suppressed and discontinuous electroencephalography activity in infants with hypoxic ischemic encephalopathy. *Neurophotonics* **3** (2016) 031408–031408.

[9] M. Darbas and S. Lohrengel, Review on mathematical modelling of electroencephalography (EEG). *Jahresber. Dtsch. Math.-Ver.* **121** (2019) 3–39.

[10] S. Pursiainen, F. Lucka and C.H. Wolters, Complete electrode model in EEG: relationship and differences to the point electrode model. *Phys. Med. Biol.* **57** (2012) 999–1017.

[11] J. Vorwerk, M. Clerc, M. Burger and C.H. Wolters, Comparison of boundary element and finite element approaches to the EEG forward problem. *Biomed. Tech.* **57** (2012) 795–798.

[12] C.H. Wolters, H. Köstler, C. Möller, J. Härdtlein, L. Grasedyck and W. Hackbusch, Numerical mathematics of the subtraction method for the modelling of a current dipole in EEG source reconstruction using finite element head models. *SIAM J. Sci. Comput.* **30** (2007) 24–45.

[13] B. Montcel, R. Chabrier and P. Poulet, Time-resolved absorption and hemoglobin concentration difference maps: a method to retrieve depth-related information on cerebral hemodynamics. *Opt Express.* **14** (2006) 12271–12287.

[14] S.R. Arridge, Optical tomography in medical imaging. *Inverse Probl.* **15** (1999) R41–R93.

[15] S.R. Arridge and J.C. Schotland, Optical tomography: forward and inverse problem. *Inverse Probl.* **25** (2009) 123010.

[16] J.P. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems. Applied Mathematical Sciences, Vol. 160. Springer (2005).

[17] A. Aubert and R. Costalat, A model of the coupling between brain electrical activity, metabolism, and hemodynamics: application to the interpretation of functional neuroimaging. *Neuroimage* **17** (2002) 1162–1181.
[18] R.B. Buxton, E.C. Wong and L.R. Frank, Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 39 (1998) 855-864.

[19] R.B. Buxton, K. Uludag, D.J. Dubowitz and T.T. Liu, Modeling the hemodynamic response to brain activation. Neuroimage 23 (2004) S220-S233.

[20] E.J. Mathias, Computational modelling of neurovascular coupling and the BOLD signal. Ph.D. Thesis, University of Canterbury, United Kingdom (2017).

[21] E.J. Mathias, A. Kenny, M.J. Plank and T. David, Integrated models of neurovascular coupling and BOLD signals: responses for varying neural activations. Neuroimage 174 (2018) 69-86.

[22] S. Sten, H. Podéus, N. Sundqvist, F. Elinder, M. Engström and G. Cedersund, A quantitative model for human neurovascular coupling with translated mechanisms from animals. PLOS Comput. Biol. 19 (2023) e1010818.

[23] F. Rapetti and G. Rousseaux, On quasi-static models hidden in Maxwell's equations. Appl. Num. Math. 79 (2014) 92-106.

[24] O. Faugeras, F. Clément, R. Deriche, R. Keriven, T. Papadopoulos, J. Roberts, T. Viéville, F. Devernay, J. Gomes, G. Hermosillo, P. Kornprobst and D. Lingrand, The inverse EEG and MEG problems: the adjoint state approach. I. The continuous case. Inria, version 1 (1999). hal.inria.fr/docs/00/07/71/12/PDF/RR-3673.pdf.

[25] IT'S Foundation, https://itis.swiss/virtual-population/tissue-properties/database/dielectric-properties/, visited on January 17th, 2025.

[26] C. Gabriel, Compilation of the Dielectric Properties of Body Tissues at RF and Microwave Frequencies (1996).

[27] H. Azizollahi, A. Aarabi and F. Wallois, Effects of uncertainty in head tissue conductivity and complexity on EEG forward modeling in neonates. Hum. Brain Mapp. 37 (2016) 3604-3622.

[28] S. Lew, D.D. Sliva, M. Choe, P.E. Grant, Y. Okada, C.H. Wolters and M.S. Hämäläinen, Effects of sutures and fontanels on MEG and EEG source analysis in a realistic infant head model. Neuroimage 76 (2013) 282-293.

[29] H. Azizollahi, M. Darbas, M.M. Diallo, A. El Badia and S. Lohrengel, EEG in neonates: forward modeling and sensitivity analysis with respect to variations of the conductivity. Math. Biosci. Eng. 15 (2018) 905-932.

[30] A. El Badia and T. Ha-Duong, An inverse source problem in potential analysis. Inverse Probl. 16 (2000) 651-663.

[31] H. Dehghani, M.E. Eames, P.K. Yalavarthy, S.C. Davis, S. Srinivasan, C.M. Carpenter, B.W. Pogue and K.D. Paulsen, Near infrared optical tomography using NIRFAST: algorithm for numerical model and image reconstruction. Commun. Numer. Methods Eng. 25 (2009) 711-732.

[32] N. Roche-Labarbe, F. Wallois, E. Ponchel, G. Kongolo and R. Grebe, Coupled oxygenation oscillation measured by NIRS and intermittent cerebral activation on EEG in premature infants. Neuroimage 36 (2007) 718-727.

[33] M. Dehaes, K. Kazemi, M. Pélégrini-Issac, R. Grebe, H. Benali and F. Wallois, Quantitative effect of the neonatal fontanel on synthetic near infrared spectroscopy measurements. Hum. Brain Mapp. 34 (2013) 878-889.

[34] A.H. Hielscher, R.E. Alcouffe and R.L. Barbour, Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues. Phys. Med. Biol. 43 (1998) 1285-1302.

[35] E. Okada and D.T. Delpy, Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer. Appl. Opt. 42 (2003) 2906-2914.

[36] S. Lohrengel, F. Mahmoudzadeh, F. Oumri, S. Salmon and F. Wallois, A homogenized cerebrospinal fluid model for diffuse optical tomography in the neonatal head. Int. J. Numer. Method Biomed. Eng. 38 (2022) e3538.

[37] S.R. Arridge and J.C. Hebden, Optical imaging in medicine. II. Modelling and reconstruction. Phys. Med. Biol. 42 (1997) 841.

[38] D. Sterratt, B. Graham, A. Gillies and D. Willshaw, Principles of Computational Modelling in Neuroscience. Cambridge University Press (2011).

[39] A.L Hodgkin and A.F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117 (1952) 500-544.

[40] R.D. Traub, J.G.R. Jefferys, R. Miles, M.A. Whittington and K. Tóth, A branching dendritic model of a rodent CA3 pyramidal neurone. J. Physiol. 481 (1994) 79-95.

[41] A. Destexhe, Z.F. Mainen and T.J. Sejnowski, Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism, J. Comput. Neurosci. 1 (1994) 195-230.

[42] H. Aurlien, I.O. Gjerde, J.H. Aarseth, G. Eldøoen, B. Karlsen, H. Skeidsvoll and N.E. Gilhus, EEG background activity described by a large computerized database. Clin. Neurophysiol. 115 (2004) 665-673.

[43] L. Kocsis, P. Herman and A. Eke, The modified Beer-Lambert law revisited. Phys. Med. Biol. 51 (2006) N91.

[44] L.F. Shampine and M.W. Reichelt, The MATLAB ODE Suite. SIAM J. Sci. Comput. 18 (1997) 1-22.
[45] F. Hecht, New development in FreeFem++ J. Numer. Math. 20 (2012) 251–266.

[46] H. Kager, W.J. Wadman and G.G. Somjen, Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations. *J. Neurophysiol.* **84** (2000) 495–512.

[47] G.B. Ermentrout and D.H. Terman, Mathematical Foundations of Neurosciences. Springer (2010).

[48] H. Jiang, Diffuse Optical Tomography. Taylor and Francis (2011).

[49] M. Nourhashemi, *Multimodal analysis of neurovascular coupling in the newborn*. Ph.D. Thesis, University of Picardie Jules Verne, France (2018).

<figure><img src="image_15.png" /></figure>

Please help to maintain this journal in open access!

This journal is currently published in open access under the [Subscribe to Open model](https://www.eds Sciences.org) (S2O). We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting [subscribers@edpsciences.org](mailto:subscribers@edpsciences.org).

More information, including a list of supporters and financial transparency reports, is available at [https://edpsciences.org/en/subscribe-to-open-s2o](https://edpsciences.org/en/subscribe-to-open-s2o).