# Mathematical modeling of photoplethysmography: model assessment and validation

Adrien Lefieux, Justine Daraize, Fabien Vergnet, Marina Vidrascu, Marie Willemet, Aniss Bendjoudi, Damiano Lombardi, Miguel Angel Fernández

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

## To cite this version:

Adrien Lefieux, Justine Daraize, Fabien Vergnet, Marina Vidrascu, Marie Willemet, et al., Mathematical modeling of photoplethysmography: model assessment and validation. 2025. hal-05134478v3

HAL Id: hal-05134478

https://inria.hal.science/hal-05134478v3

Preprint submitted on 29 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mathematical modeling of

photoplethysmography: model assessment and

validation

Adrien Lefieux¹, Justine Daraize², Fabien Vergnet¹,
Marina Vidrascu¹, Marie Willemet², Aniss Bendjoudi²,
Damiano Lombardi¹, Miguel A. Fernández¹*

¹Sorbonne Université, CNRS, LJLL, Inria, Paris, France.

²Withings SAS, Issy-les-Moulineaux, France.

*Corresponding author(s). E-mail(s): miguel.fernandez@inria.fr;
Contributing authors: adrien.lefieux@inria.fr;
justine.daraize@withings.com; fabien.vergnet@sorbonne-universite.fr;
marina.vidrascu@inria.fr; marie.willemet@withings.com;
aniss.bendjoudi@withings.com; damiano.lombardi@inria.fr;

**Abstract**

**Purpose.** Photoplethysmography (PPG) is a well-known technique employed to assess optically perfused bio-tissue volume changes. A PPG apparatus consists of a light emitter and a receptor. The analysis of the received light is used to infer properties of the illuminated tissues. This paper aims at presenting a novel distributed mathematical model for PPG signals, which combines a poroelastic model of tissue perfusion with a diffusion model for light absorption and scattering.

**Methods.** We assume that the tissues undergo small deformations, allowing for a linear poroelastic description of perfusion. Since many PPG devices are applied to the fingertips (due to the rich vascularization in that area) we model the system specifically on the finger, with arterial blood pulse pressure serving as the primary perfusion driver. The numerical discretization of the governing equations is carried out using the finite element method. After calibration of the model, 216k simulations are performed with varying parameters (quasi-Monte
Carlo approach). The aggregated results for two key biomarkers, AC/DC PPG amplitude ratio and pulse pressure, are compared against experimental PPG and pulse pressure measurements obtained from 20 volunteers.

**Results.** Within the prescribed parameter ranges, numerical simulations successfully reproduce the AC/DC PPG amplitude biomarker in both the red and infrared wavelengths, with a few outliers observed for the green. Furthermore, in the vicinity of the combined red and infrared measured biomarkers, there are simulated biomarkers whose corresponding pulse pressures closely match the mapped-to-finger measured pulse pressure, with a difference of less than 1 mmHg.

**Conclusion.** This work demonstrates the relevance of the proposed mathematical model for simulating PPG signals and highlights its potential for estimating tissue perfusion parameters, such as arterial pulse pressure.

**Keywords:** Photoplethysmography (PPG) signal, biomarker, poroelasticity, diffusion bio-optics, mathematical modeling, numerical simulation, finite element method.

# 1 Introduction

Photoplethysmography (PPG) is a non-invasive technique that relies on light propagation for blood volume change assessment, in particular through tissues in the visible and near infrared wavelengths. PPG is widely used in applications such as heart rate monitoring (commonly found in smartwatches) and for measuring arterial blood oxygen saturation, typically using fingertip clips (see, e.g., [2, 34] and the references therein).

A PPG apparatus is essentially a light emitter and a photo-receptor, with the received light used to infer properties (e.g., mechanical or physiological) of the illuminated tissues or the system subtending to it. While the precise origin of the PPG signal is still debated (see, e.g., [35]), it is commonly believed that the pulsatility of the PPG signal arises, in part, from mechanical changes in blood volume and deformations of perfused tissues caused by cardiac activity, particularly due to blood pressure variations.

Understanding the relationship between inputs and outputs in the PPG apparatus-tissue system is key to effectively interpreting and utilizing PPG signals. There
is now a large body of research on the computer simulation of photoplethysmography
signals (see, e.g., [35] and references therein). Nevertheless, as regards the relationship
between blood pressure and the PPG signal, much ongoing research involves machine
learning techniques which, most often, lack proper physiological interpretability (see,
e.g., [13, 19]). Alternatively, common physiological models rely on simple laws, such
as the Beer-Lambert law for light propagation, or the zero-dimensional Windkessel
model for arterial networks (see, e.g., [38]). To the best of our knowledge, physiologi-
cally meaningful combined opto-mechanical models for PPG are scarce. In this work,
we propose to advance PPG modeling by introducing and evaluating a distributed
mathematical model that combines a diffusion approximation for light absorption and
scattering with a poroelastic model of tissue perfusion.

The Beer-Lambert law is particularly applicable in cases of strong absorption
and weak scattering. However, in the visible and near-infrared wavelengths, biological
tissues exhibit strong scattering relative to absorption (see, e.g., [53]). Fortunately,
in such situations, particularly in the case of skin tissue (see [51]), it is justifiable to
employ the so-called diffusion approximation of the radiative transfer equation (RTE)
(see [29, 53]). This hence allows for solving a 3D quasi-static optical diffusion problem,
instead of the 5D RTE. With such an approach, the absorption parameter is spatially
dependent and, based on assumptions regarding the origin of the PPG signal, it is
modulated by the pulsatile blood volume fraction in the perfused tissue. As regards
blood tissue perfusion, we choose to model the illuminated medium as a homogenized
poroelastic material (see, e.g., [12, 15]). More precisely, we assume that the perfused
tissues undergo small deformations, allowing the use of a first order approximation of
the deformation and a linear poroelastic model (see, e.g., [4, 5, 11]). Note that such a
model accounts for both fluid and solid inertial effects as well as their incompressibility.
Importantly, the perfusion arterial pulse pressure is accounted here at the boundary
of the domain. A preliminary evaluation of the proposed PPG mathematical model
is conducted by simulating a range of parameter values and comparing the results to
PPG data collected from a cohort of 20 volunteers (see [16]).

This paper is organized as follows. Section 2 presents first the poroelastic and opti-
cal models separately, followed by the formulation of the combined poromechanical-
optical mathematical model for PPG. Details on the numerical discretization of the
resulting systems are also provided. The section is concluded with a description of
the experimental PPG dataset used for the assessment and validation of the proposed
PPG simulator. Section 3 begins with a qualitatively assessment of the model after
calibration, via comparisons against the PPG measurements on a single heart beat.
Next, a broader range of model parameters is explored to evaluate the simulations'
ability to reproduce the entire PPG dataset, focusing on a specific biomarker. Finally,
the results, limitations and perspectives of this work are discussed in Section 4.

## 2 Materials and Methods

In this section, we first present the fundamental ingredients of PPG measurements.
Then, on such bases, we propose an opto-poromechanical model of PPG signals. As
already discussed in the introduction, a PPG system is composed of two elements: a
light emitting apparatus, such as a LED, and a light receptor, such as a photodiode.
The dominant hypothesis for the origin of the photoplethysmographic signal, as its
name suggests, is that a variation in blood volume content in the illuminated perfused
tissue induces a change in light flux received by the photoreceptor, as photons are
absorbed, primarily, by hemoglobin.

The considered geometrical description of the system under study is presented in
Figure 1, which corresponds to a reflectance PPG measurement on the finger (i.e.,
the light source and photodetector are on the same side of the skin). Nevertheless,
the proposed model can be extended to other body locations or different acquisition
settings. The finger domain $\Omega \subset \mathbb{R}^3$ is assumed to be made of two components $\Omega_T$
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Schematic description of a PPG optoelectronics apparatus by reflectance applied on the phalanx region $\Omega \stackrel{\text{def}}{=} \Omega_T \cup \Omega_B \cup \Gamma_B$ where, $\Omega_T$ and $\Omega_B$ respectively correspond to the tissue and bone regions of the finger and $\Gamma_B$ denotes the tissue-bone interface. The lateral side $\Gamma_T^e \cup \Gamma_B^e$ represents the (artificial) external boundary of the computational domain and $\Gamma_S$ the skin boundary. Finally, $\Gamma_E$, $\Gamma_R$ and $\Gamma_D$ respectively denote the skin zones in contact with the emitter, the receptor and the remainder part of the PPG device. Emitted photons through $\Gamma_E$ are absorbed by the perfused tissue in $\Omega$, primarily hemoglobin, and the variation in exiting light flux on $\Gamma_R$ is recorded by the photoreceptor.</figcaption></figure>

and $\Omega_B$, which respectively denote the regions filled by the finger tissue and the bone. Furthermore, we assume that the boundary of the bone domain is partitioned as $\partial\Omega_B = \Gamma_B \cup \Gamma_B^e$, where $\Gamma_B$ represents the tissue-bone interface $\Gamma_B \stackrel{\text{def}}{=} \partial\Omega_T \cap \partial\Omega_B$, and $\Gamma_B^e$ stands for the (artificial) external bone boundary. Similarly, the tissue boundary is partitioned as $\partial\Omega_T = \Gamma_S \cup \Gamma_T^e \cup \Gamma_B \cup \Gamma_E \cup \Gamma_R \cup \Gamma_D$. Here, $\Gamma_S$ corresponds to the skin boundary, $\Gamma_T^e$ to the tissue (artificial) external boundary, while $\Gamma_E$, $\Gamma_R$ and $\Gamma_D$ respectively denote the skin zones in contact with the emitter, the receptor and the remainder device boundary.

Accordingly, the proposed model is composed of two essential components: a mechanical (poroelastic) model for the perfused tissues' deformation under blood pressure, combined with an optical model for light propagation within the tissues (diffusion), whose absorption parameter depends on the local blood volume content.
## 2.1 Poroelastic Model

We consider here the linear poroelastic model reported in [4, 5], obtained by linearization of the fully non-linear coupled poroelastic system introduced in [12]. The model generalizes the so-called Biot theory of consolidation (see, e.g., [6, 15]), by accounting for both solid and fluid inertia (critical when the fluid and the solid have comparable densities, as in biological tissues) and viscous effects in the fluid phase (possibly in the solid one as well).

More precisely, we assume that the finger tissue domain $\Omega_T$ is filled by an incompressible porous medium involving an homogeneous elastic matrix, of volume fraction $1 - \phi_0$, $\phi_0 \in (0, 1)$, and an homogeneous, viscous and Newtonian fluid, of volume fraction $\phi_0$. The resulting model (in the absence of body forces) reads as follows: Find the solid displacement $\mathbf{d}: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$, the solid velocity $\dot{\mathbf{d}}: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$, the fluid velocity $\mathbf{u}: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$ and the tissue interstitial pressure $p: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}$ with $\dot{\mathbf{d}} = \partial_t \mathbf{d}$ and such that

$$ \left\{ \begin{array}{ll} \rho_s(1 - \phi_0)\partial_t \dot{\mathbf{d}} - \mathbf{div}\sigma_s(\mathbf{d}) + \frac{\phi_0^2}{k}(\dot{\mathbf{d}} - \mathbf{u}) + (1 - \phi_0)\nabla p = \mathbf{0} & \text{in } \Omega_T, \\ \\ \rho_f\phi_0\partial_t \mathbf{u} - \mathbf{div}(\phi_0\sigma_f(\mathbf{u})) + \frac{\phi_0^2}{k}(\mathbf{u} - \dot{\mathbf{d}}) + \phi_0\nabla p = \mathbf{0} & \text{in } \Omega_T, \\ \\ \mathbf{div}((1 - \phi_0)\dot{\mathbf{d}} + \phi_0\mathbf{u}) = 0 & \text{in } \Omega_T. \end{array} \right. \quad (1) $$

Here, the first and second relations respectively express the momentum conservation of the elastic and fluid phases, while the third enforces total mass conservation. The symbols $\rho_s$ and $\rho_f$ respectively denote the solid and fluid densities [g/cm³], while the Cauchy stress tensors $\sigma_s$ and $\sigma_f$ are given by the relations

$$ \sigma_i(\mathbf{g}) \stackrel{\text{def}}{=} \lambda_i(\mathbf{div}\mathbf{g})\mathbf{I} + 2\mu_i\epsilon(\mathbf{g}) \quad i \in \{\text{s, f}\}, \quad \epsilon(\mathbf{g}) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla\mathbf{g} + (\nabla\mathbf{g})^\text{T}), $$
where $I$ is the identity tensor, $\lambda_s$ and $\mu_s$ [Ba] are the two Lamé coefficients of the solid phase and $\mu_f$ and $\lambda_f$ are, respectively, the first (shear) and second dynamic viscosity coefficients of the fluid [Bas]. Finally, the coefficient $k$ is the (isotropic) hydraulic conductivity [cm$^2$ Ba$^{-1}$ s$^{-1}$].

The poroelastic system (1) has to be complemented with appropriate initial and boundary conditions. For the sake of simplicity, we assume that the system is initially at rest, i.e., $d(\cdot, 0) = 0$, $\dot{d}(\cdot, 0) = 0$ and $u(\cdot, 0) = 0$ in $\Omega_T$. As regards the boundary conditions, it is useful to introduce the total solid and fluid tractions, $t_s$ and $t_f$ respectively, on $\partial\Omega_T$ (see [11]):

$$ \begin{cases} t_s(d, p) \stackrel{\text{def}}{=} \sigma_s(d)n - (1 - \phi_0)pn, \\ t_f(u, p) \stackrel{\text{def}}{=} \phi_0(\sigma_f(u)n - pn), \end{cases} $$

where $n$ stands for the outward normal on $\partial\Omega_T$. In what follows, we let $\{\tau_1, \tau_2\}$ denote the local tangent basis on $\partial\Omega_T$. The following boundary conditions are enforced on the solid and fluid phases:

$$ \left\{ \begin{array}{l} d = 0 \\ u = 0 \end{array} \right. \quad \text{on} \quad \Gamma_B \cup \Gamma_D \cup \Gamma_E \cup \Gamma_R, \qquad (2) $$

$$ \left\{ \begin{array}{l} d \cdot n = 0 \\ t_s(d, p) \cdot \tau_i = 0 \\ t_f(u, p) = -\phi_0 P n \end{array} \right. \quad \text{on} \quad \Gamma_T^e, \qquad (3) $$

$$ \left\{ \begin{array}{l} u = \dot{d} \\ t_s(d, p) + t_f(u, p) = 0 \end{array} \right. \quad \text{on} \quad \Gamma_S, \qquad (4) $$

with $i = 1, 2$. The Dirichlet conditions (2), on the boundary in contact with the bone and the PPG device $\Gamma_B \cup \Gamma_D \cup \Gamma_E \cup \Gamma_R$, respectively enforce that the solid phase is fixed
and that the fluid phase cannot slip. On the artificial inlet boundary $\Gamma_T^e$, a symmetric condition is prescribed on the solid phase, while, for the fluid, we enforce a normal traction in terms of a given inlet temporal pressure profile $P: \mathbb{R}^+ \to \mathbb{R}$. Finally, on the (remaining) skin boundary $\Gamma_S$, the relations (4) impose a total kinematic and dynamic balance between the fluid and solid phases.

The blood volume fraction in the deformed configuration $\phi: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}$ can be reconstructed in terms of the solid displacement $d$ (at first-order, see [5, Appendix A]) as

$$ \phi \approx \phi_0 + \operatorname{div} d, $$

so that the pulsatile blood volume fraction $\delta\phi$ is given by

$$ \delta\phi \stackrel{\text{def}}{=} \operatorname{div} d. \qquad (5) $$

## 2.2 Bio-optics Model

We employ the so-called diffusion approximation of the radiative transfer equation (RTE) for light scattering and absorption (see, e.g., [29, 53]). Before presenting the complete optical mathematical model, we first introduce some basic notations and several essential radiometric quantities. In what follows, the symbol $\mathbb{S}^2$ denotes the unit sphere in $\mathbb{R}^3$. The *radiance* field $L: \Omega \times \mathbb{S}^2 \times \mathbb{R}^+ \to \mathbb{R}$ provides the amount of light power $L(x, s, t)$ (per unit area and per unit solid angle [$mWm^{-2}sr^{-1}$]) flowing at a point $x \in \Omega$, in a specific direction $s \in \mathbb{S}^2$ and at a given time instant $t \ge 0$. For the purpose of the present work, the fundamental radiometric quantity is the *fluence* rate $\Phi: \Omega \times \mathbb{R}^+ \to \mathbb{R}$, which represents the total light power per unit area [$mW/cm^2$] arriving at a point $x \in \Omega$ from all directions in $\mathbb{S}^2$, viz.,

$$ \Phi(x, t) \stackrel{\text{def}}{=} \int_{\mathbb{S}^2} L(x, s, t)ds. $$
<figure><img src="image_4.png" /><figcaption>Fig. 2 Half-spheres separated by the boundary $\Sigma$. The normal $n$ is oriented outward from the tissue.</figcaption></figure>

The diffusion approximation of the RTE is derived by combining a first-order spherical harmonics decomposition of the radiance $L$ with Fick's law, which expresses the first moment of the radiance $L$ over $s$ as a linear function of the gradient of the fluence rate $\Phi$ (see, e.g., [29, 53]), namely,

$$L(x, s, t) \approx \frac{1}{4\pi}\Phi(x, t) - \frac{3}{4\pi}D\nabla\Phi(x, t) \cdot s, \quad (6)$$

with the diffusion coefficient $D$ given by

$$D \stackrel{\text{def}}{=} \frac{1}{3(\mu_a + \mu'_s)}. \quad (7)$$

Here, the symbol $\mu_a$ denotes the absorption coefficient [cm$^{-1}$] and $\mu'_s$ is the reduced scattering coefficient [cm$^{-1}$]. Furthermore, owing to the much faster speed of light with respect to the pulse wave, we consider the quasi-static diffuse approximation equation for light absorption and scattering in the finger phalanx domain $\Omega$ (see Figure 1): Find the fluence rate $\Phi: \Omega \times \mathbb{R}^+ \to \mathbb{R}$ such that

$$-\mathrm{div}(D\nabla\Phi) + \mu_a\Phi = 0 \quad \text{in } \Omega. \quad (8)$$
In this paper, we will consider a simple monolayer skin model in which the reduced scattering coefficient $\mu'_s$ is assumed to be piece-wise constant in $\Omega$, that is,

$$ \mu'_s = \begin{cases} \mu'_{s,T} & \text{in } \Omega_T, \\ \mu'_{s,B} & \text{in } \Omega_B, \end{cases} $$

with the values $\mu'_{s,T}$ and $\mu'_{s,B}$ depending only on the light wavelength $\lambda$, typically through power-laws of the form $\mu'_s(\lambda) = \alpha\lambda^\beta$ (accounting for Mie and Rayleigh scattering see, e.g., [30, 31]). In the bone region $\Omega_B$, the absorption coefficient $\mu_a$ will also be assumed to be constant and function only of the wave length $\mu_a = \mu_{a,B}(\lambda)$ (see, e.g., [49]). In the tissue region $\Omega_T$ however, the coefficient $\mu_a$ will depend on both space and time, as a function of the pulsatile blood volume content $\delta\phi$, given by (5). This point will be further developed in Section 2.3.1 below.

**Remark 1.** Recently, several studies have argued that $D \approx (3\mu'_s)^{-1}$ (see, e.g., [18] and the references therein). Such a point is not considered here. Additionally, according to [32], a rule of thumb for the validity of the diffusion approximation of the RTE equation is that $\mu'_s/\mu_a > 10$. Furthermore, here, we assume that the refractive index is homogeneous within $\Omega$.

While boundary conditions for elliptic problems such as (8) are standard (see, e.g., [21]), ensuring that they correctly approximate the incoming radiance boundary condition of the RTE requires some care (see, e.g., [25, 45]). Let $\Sigma \subset \partial\Omega$ be an arbitrary part of $\partial\Omega$ with outward pointing normal $n$. At each point of $\Sigma$, we define the inward and outward unit hemispheres $S_-^2 \stackrel{\text{def}}{=} \{s \in S^2 : s \cdot n < 0\}$ and $S_+^2 \stackrel{\text{def}}{=} \{s \in S^2 : s \cdot n > 0\}$, respectively (see Figure 2). Boundary conditions for (8) are generally
derived from the following incoming flux balance condition:

$$ \int_{S_{\perp}} L(\cdot, \mathbf{s}, t) \mathbf{s} \cdot \mathbf{n} \, d\mathbf{s} = R \int_{S_{\perp}} L(\cdot, \mathbf{s}, t) \mathbf{s} \cdot \mathbf{n} \, d\mathbf{s} + l \quad \text{on } \Sigma, \qquad (9) $$

where R denotes the effective reflection coefficient and $l: \Sigma \times \mathbb{R}^+ \to \mathbb{R}$ stands for an external light source (in $[\text{mW/cm}^2]$). The relation (9) enforces that the total inward flux equals the sum of the reflected part of the outward flux and the external incoming flux. By inserting the diffusion approximation (6) into (9), we get

$$ \frac{\Phi}{4} + \frac{D}{2} \frac{\partial \Phi}{\partial \mathbf{n}} = R \left( \frac{\Phi}{4} - \frac{D}{2} \frac{\partial \Phi}{\partial \mathbf{n}} \right) + l \quad \text{on } \Sigma, $$

which yields the following general Robin boundary condition for (8):

$$ D \frac{\partial \Phi}{\partial \mathbf{n}} + \frac{\Phi}{2C_R} = \frac{2l}{1+R} \quad \text{on } \Sigma, \qquad (10) $$

where

$$ C_R \stackrel{\text{def}}{=} \frac{1+R}{1-R}. $$

On the parts of the boundary which are not in contact with air, viz., $\Gamma_T^e \cup \Gamma_B^e \cup \Gamma_D \cup \Gamma_E \cup \Gamma_R$ (see Figure 1), we will assume that incident light is fully transmitted, so that $R = 0$ and $C_R = 1$. On the remainder boundary part $\Gamma_S$ (i.e., with air as external medium), an approximation of R is provided by the following function (see, e.g., [25]): $R(n_-) \approx -1.4399n_-^2 + 0.7099n_-^1 + 0.6681 + 0.0636n_-$, where $n_-$ denotes the internal refraction index. In our application, $n_- \approx 1.4$ (biological tissue with 65% water content, see [31]), which yields $R \approx 0.5$ (around half of the incident light is reflected into the tissue) and $C_R \approx 3$. Therefore, owing to (10), the light diffusion
equation (8) is complemented with the following boundary conditions on $\partial\Omega$:

$$ \left\{ \begin{array}{ll} D \frac{\partial \Phi}{\partial n} + \frac{\Phi}{6} = 0 & \text{on } \Gamma_S, \\ D \frac{\partial \Phi}{\partial n} + \frac{\Phi}{2} = 0 & \text{on } \Gamma_T^e \cup \Gamma_B^e \cup \Gamma_R \cup \Gamma_D, \\ D \frac{\partial \Phi}{\partial n} + \frac{\Phi}{2} = 2l & \text{on } \Gamma_E, \end{array} \right. \quad (11) $$

where $l$ denotes the irradiant flux of the emitter.

## 2.3 Combined Poromechanical-Optical PPG Model

PPG signals are generated by variations in tissue perfusion, which modulate the amount of light reaching the photodetector through hemoglobin absorption. Therefore, the remaining modeling ingredient, in the poroelastic and optical systems introduced above, is the description of how changes in the pulsatile blood volume fraction $\delta\phi$ induce variations of the absorption coefficient $\mu_a$ in the tissue region $\Omega_T$.

### 2.3.1 Blood Flow Induced Changes in Light Absorption

Following [30, 35, 44], we describe the absorption coefficient $\mu_a$ in the skin tissue region $\Omega_T$ as the sum of the distinct chromophores (light-absorbing molecules) that compose the skin layers:

$$ \mu_a \stackrel{\text{def}}{=} (S\phi_0 + \delta\phi)\mu_{a,H_bO_2} + (1-S)\phi_0\mu_{a,H_b} + (1-(\phi_0 + \delta\phi))\mu_{a,T} \quad \text{in } \Omega_T, \quad (12) $$

where $\mu_{a,H_bO_2}$, $\mu_{a,H_b}$ and $\mu_{a,T}$ denote the absorption coefficient of the oxygenized hemoglobin, deoxygenized hemoglobin and bloodless skin tissues, respectively, and $S$ stands for the arterio-venous blood oxygen. Note that, in (12), absorption variations are not induced by changes in deoxygenized hemoglobin. Furthermore, the parameters
$\mu_{a,H_b O_2}$ and $\mu_{a,H_b}$ depend on the hematocrit $H$, via the relation

$$ \mu_{a,i} \stackrel{\text{def}}{=} \ln(10) H c_{Hb} \epsilon_{a,i}, \quad i \in \{\text{H}_b\text{O}_2, \text{Hb}\}, \qquad (13) $$

where $\epsilon_{a,i}$ denotes the extinction coefficient ($[\text{L cm}^{-1} \text{mol}^{-1}]$) and $c_{Hb} = 5.1713 \cdot 10^{-3} \text{mol/L}$ the hemoglobin molar concentration (see [43, 50]). Note that $\epsilon_{a,H_b O_2}$, $\epsilon_{a,H_b}$ (and so $\mu_{a,H_b O_2}$, $\mu_{a,H_b}$) and $\mu_{a,T}$ will be functions of the light wave length $\lambda$.

### 2.3.2 PPG Signal Mathematical Model

By collecting all the above modeling ingredients, we can now summarize the main components of the PPG mathematical model introduced in this paper.

#### Skin Tissue Perfusion Model

Owing to (1)–(4), the first fundamental ingredient of the proposed PPG model is the poroelastic description of the tissue perfusion in the skin domain $\Omega_T$. Let $P: \mathbb{R}^+ \to \mathbb{R}$ be a given temporal perfusion pressure profile on $\Gamma_T^e$ and assume that all the bio-physical parameters reported in Table 1 are known. We proceed by solving the

<table><caption>**Table 1** Bio-physical parameters of the poroelastic model.</caption><thead><tr><th>Domain</th><th colspan="8">Ω<sub>T</sub></th></tr><tr><th>Parameters</th><th>φ<sub>0</sub></th><th>k</th><th>ρ<sub>s</sub></th><th>λ<sub>s</sub></th><th>μ<sub>s</sub></th><th>ρ<sub>f</sub></th><th>λ<sub>f</sub></th><th>μ<sub>f</sub></th></tr><tr><th>Units</th><th>%</th><th>cm<sup>2</sup> Ba<sup>-1</sup> s<sup>-1</sup></th><th>g cm<sup>-3</sup></th><th>Ba</th><th>Ba</th><th>g cm<sup>-3</sup></th><th>Ba s</th><th>Ba s</th></tr></thead></table>

following problem: Find the solid phase displacement $d: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$, the solid phase velocity $\dot{d}: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$, the fluid phase velocity $u: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}^3$ and the
interstitial pressure $p: \Omega_T \times \mathbb{R}^+ \to \mathbb{R}$, with $\dot{d} = \partial_t d$ and such that

$$
\left\{
\begin{array}{ll}
\rho_s(1 - \phi_0)\partial_t \dot{d} - \mathbf{div}\sigma_s(d) + \frac{\phi_0^2}{k}(\dot{d} - u) + (1 - \phi_0)\nabla p = \mathbf{0} & \text{in } \Omega_T, \\
\\
\rho_f\phi_0\partial_t u - \mathbf{div}(\phi_0\sigma_f(u)) + \frac{\phi_0^2}{k}(u - \dot{d}) + \phi_0\nabla p = \mathbf{0} & \text{in } \Omega_T, \\
\\
\mathbf{div}((1 - \phi_0)\dot{d} + \phi_0 u) = 0 & \text{in } \Omega_T, \\
\\
d = \mathbf{0}, \quad u = \mathbf{0} & \text{on } \Gamma_D \cup \Gamma_E \cup \Gamma_R \cup \Gamma_B, \\
\\
d \cdot n = 0, \quad t_s(d, p) \cdot \tau_i = 0, \quad t_f(u, p) = -\phi_0 P n & \text{on } \Gamma_T^e, \\
\\
\dot{d} = u, \quad t_s(d, p) + t_f(u, p) = \mathbf{0} & \text{on } \Gamma_S, \\
\\
d(\cdot, 0) = \mathbf{0}, \quad \dot{d}(\cdot, 0) = \mathbf{0}, \quad u(\cdot, 0) = \mathbf{0} & \text{in } \Omega_T.
\end{array}
\right.
\tag{14}
$$

The pulsatile blood volume fraction $\delta\phi : \Omega_T \times \mathbb{R}^+ \to \mathbb{R}$ can thereafter be directly recovered from the displacement field $d$ via the relation (5), viz., $\delta\phi = \text{div } d$.

Optical Model

The second key component of the PPG model consists in the description of the light
propagation within the finger phalanx domain Ω, via the relations (8), (11) and (12).
Let l : R+ → R be a given irradiant flux on the emitter ΓE and assume that all
the bio-physical parameters reported in Table 2 are known, the reduced scattering,

<table><caption>**Table 2** Bio-physical parameters of the optical model.</caption>
  <thead>
    <tr>
      <th>Domains</th>
      <th colspan="6">Ω<sub>T</sub></th>
      <th colspan="2">Ω<sub>B</sub></th>
    </tr>
    <tr>
      <th>Parameters</th>
      <th>ε<sub>a,H<sub>b</sub>O<sub>2</sub></sub></th>
      <th>ε<sub>a,H<sub>b</sub></sub></th>
      <th>μ<sub>a,T</sub></th>
      <th>μ'<sub>s,T</sub></th>
      <th>S</th>
      <th>H</th>
      <th>μ<sub>a,B</sub></th>
      <th>μ'<sub>s,B</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Units</th>
      <td>L cm<sup>-1</sup> mol<sup>-1</sup></td>
      <td>L cm<sup>-1</sup> mol<sup>-1</sup></td>
      <td>cm<sup>-1</sup></td>
      <td>cm<sup>-1</sup></td>
      <td>%</td>
      <td>%</td>
      <td>cm<sup>-1</sup></td>
      <td>cm<sup>-1</sup></td>
    </tr>
  </tbody>
</table>

absorption and diffusion parameters coefficients can be evaluated, from (13) and (14),
via the relations:

$$
\begin{align}
\mu_a &= \begin{cases}
(S\phi_0 + \operatorname{div} \boldsymbol{d}) \mu_{a,H_b O_2} + (1-S)\phi_0\mu_{a,H_b} + (1-(\phi_0 + \operatorname{div} \boldsymbol{d}))\mu_{a,T} & \text{in } \Omega_T, \\
\mu_{a,B} & \text{in } \Omega_B,
\end{cases} \nonumber \\
\mu'_s &= \begin{cases}
\mu'_{s,T} & \text{in } \Omega_T, \\
\mu'_{s,B} & \text{in } \Omega_B,
\end{cases} \nonumber \\
D &= \frac{1}{3(\mu_a + \mu'_s)} \quad \text{in } \Omega. \tag{15}
\end{align}
$$

We then proceed by solving the following quasi-static light diffusion problem: Find
the fluence rate $\Phi : \Omega \times \mathbb{R}^+ \rightarrow \mathbb{R}$ such that

$$
\left\{
\begin{array}{ll}
-\operatorname{div}(D\nabla\Phi) + \mu_a\Phi = 0 & \text{in } \Omega \\
D\frac{\partial\Phi}{\partial n} + \frac{\Phi}{6} = 0 & \text{on } \Gamma_S, \\
D\frac{\partial\Phi}{\partial n} + \frac{\Phi}{2} = 0 & \text{on } \Gamma_T^e \cup \Gamma_B^e \cup \Gamma_R \cup \Gamma_D, \\
D\frac{\partial\Phi}{\partial n} + \frac{\Phi}{2} = 2l & \text{on } \Gamma_E.
\end{array}
\right.
\tag{16}
$$

PPG Signal Generation

Finally, the PPG signal recorded by the photodiode (Lambertian) corresponds to the
total outgoing light flux $Q_R : \mathbb{R}^+ \rightarrow \mathbb{R}$ on $\Gamma_R$, which can simply be evaluated from $\Phi$
via the following expression (assuming similar refractive indices on both sides of $\Gamma_R$,
see Section 2.2):

$$
\int_{\Gamma_R} \int_{S_+^2} Ls \cdot n \, ds \approx Q_R \stackrel{\text{def}}{=} \int_{\Gamma_R} \left( \frac{\Phi}{4} - \frac{D}{2} \frac{\partial \Phi}{\partial n} \right). \quad (17)
$$
## 2.4 Numerical Method for PPG Simulation

This section describes the computational methods employed to numerically solve the continuous problems (14) and (16), which thereafter enable the simulation of the PPG signal tanks to (17).

The discretization in time is performed with a finite difference scheme, while a finite element approach is used for the spatial discretization. In what follows, we consider a partition of the time interval $(0, T)$ made of $N \in \mathbb{N}^*$ sub-intervals $\{(t_{n-1}, t_n)\}_{1 \le n \le N}$, where $t_n \stackrel{\text{def}}{=} \tau n$, $T = t_N$, and $\tau > 0$ denotes the time-step length. We shall also make use of the following standard notation for the first-order backward finite difference:

$$ \partial_{\tau} f^n \stackrel{\text{def}}{=} \frac{1}{\tau} (f^n - f^{n-1}). $$

In order to ease the presentation, we assume that $\Omega_T$ and $\Omega_B$ are polyhedral domains and that the artificial exterior boundary $\Gamma_T^e \cup \Gamma_B^e$ is flat. We let $\mathcal{T}_h^{\Omega_T}$ and $\mathcal{T}_h^{\Omega_B}$ be given triangulations of $\Omega_T$ and $\Omega_B$, respectively, where the subindex $h > 0$ denotes the level of mesh refinement. We will only consider triangulations $\mathcal{T}_h^{\Omega_T}$ and $\mathcal{T}_h^{\Omega_B}$ that are fitted on $\Gamma_B$, so that $\mathcal{T}_h \stackrel{\text{def}}{=} \mathcal{T}_h^{\Omega_T} \cup \mathcal{T}_h^{\Omega_B}$ yields a triangulation of the whole domain $\Omega$.

### 2.4.1 Poroelastic Model Discretization

The poroelastic system (14) is discretized in time using a backward-Euler scheme. Piece-wise affine continuous finite elements are used for the spatial discretization. We hence consider the following finite element spaces for the fluid and solid phases

$$ \mathbf{V}_h^f \stackrel{\text{def}}{=} \{ \boldsymbol{v}_h \in [C^0(\overline{\Omega_T})]^3 : \boldsymbol{v}_h|_K \in [\mathbb{P}_1]^3 \quad \forall K \in \mathcal{T}_h^{\Omega_T}, \quad \boldsymbol{v}_h|_{\Gamma_D \cup \Gamma_E \cup \Gamma_R \cup \Gamma_B} = \mathbf{0} \}, $$

$$ \mathbf{V}_h^s \stackrel{\text{def}}{=} \{ \boldsymbol{v}_h \in \mathbf{V}_h^f : \boldsymbol{v}_h \cdot \boldsymbol{n}|_{\Gamma_T^e} = 0 \}, $$
and for the interstitial pressure

$$Q_h \stackrel{\text{def}}{=} \{ q_h \in C^0(\overline{\Omega_T}) : q_h|_{K} \in \mathbb{P}_1, \forall K \in \mathcal{T}_h^{\Omega_T} \}.$$

Alternative spatial approximations have been advocated in [4, 5], by building on standard inf-sup stable velocity/pressure pairs for the Stokes problem. Here, for simplicity, we consider equal-order affine discretizations in all the unknowns. In order to cope with the lack of inf-sup stability of the triplet $V_h^s \times V_h^f \times Q_h$ (see, e.g., [4]), we consider the following weakly consistent pressure-stabilization term, from the Brezzi-Pitkäranta method for the Stokes equations (see [8]):

$$s_h(p_h, q_h) \stackrel{\text{def}}{=} \frac{\gamma_p h^2}{\mu_f} \int_{\Omega} \nabla p_h^n \cdot \nabla q_h \quad \forall p_h, q_h \in Q_h, \quad (18)$$

where $\gamma_p > 0$ is a dimensionless user-defined parameter. This guarantees the discrete inf-sup stability condition (A1) reported in Appendix A. Additionally, a Nitsche's type mortaring is utilized to weakly enforce the kinematic and dynamic coupling between the fluid and solid phases on $\Gamma_S$ (see, e.g., [10, 11, 26]).

To wrap up, the proposed fully discrete approximation of (14) reads as follows: For $1 \le n \le N$, find $(\mathbf{d}_h^n, \dot{\mathbf{d}}_h^n, \mathbf{u}_h^n, p_h^n) \in V_h^s \times V_h^s \times V_h^f \times Q_h$ with $\dot{\mathbf{d}}_h^n = \partial_\tau \mathbf{d}_h^n$ and such that

$$\begin{aligned}
& \int_{\Omega_T} \rho_s(1 - \phi_0) \partial_\tau \dot{\mathbf{d}}_h^n \cdot \mathbf{v}_h^s + \int_{\Omega_T} \boldsymbol{\sigma}_s(\mathbf{d}_h^n) \cdot \boldsymbol{\varepsilon} (\mathbf{v}_h^s) + \int_{\Omega_T} \frac{\phi_0^2}{k} (\dot{\mathbf{d}}_h^n - \mathbf{u}_h^n) \cdot (\mathbf{v}_h^s - \mathbf{v}_h^f) \\
& \quad + \int_{\Omega_T} \rho_f \phi_0 \partial_\tau \mathbf{u}_h^n \cdot \mathbf{v}_h^f + \int_{\Omega_h} \phi_0 \boldsymbol{\sigma}_f(\mathbf{u}_h^n) \cdot \boldsymbol{\varepsilon}(\mathbf{v}_h^f) + s_h(p_h^n, q_h) \\
& \quad - \int_{\Omega_T} p_h^n \operatorname{div} \left( (1 - \phi_0) \mathbf{v}_h^s + \phi_0 \mathbf{v}_h^f \right) + \int_{\Omega_T} q_h \operatorname{div} \left( (1 - \phi_0) \dot{\mathbf{d}}_h^n + \phi_0 \mathbf{u}_h^n \right) \quad (19) \\
& \quad - \int_{\Gamma_S} \mathbf{t}_f(\mathbf{u}_h^n, p_h^n) \cdot (\mathbf{v}_h^f - \mathbf{v}_h^s) - \int_{\Gamma_S} q_h(\mathbf{u}_h^n - \dot{\mathbf{d}}_h^n) \cdot \mathbf{n} \\
& \quad + \frac{\gamma_n \phi_0 \mu_f}{h} \int_{\Gamma_S} (\mathbf{u}_h^n - \dot{\mathbf{d}}_h^n) \cdot (\mathbf{v}_h^f - \mathbf{v}_h^s) = - \int_{\Gamma_T} \phi_0 P(t_n) \mathbf{n} \cdot \mathbf{v}_h^f
\end{aligned}$$
for all $(v_h^s, v_h^f, q_h) \in V_h^s \times V_h^f \times Q_h$ and where the symbol $\gamma_n > 0$ denotes the so-called Nitsche penalty parameter (dimensionless).

### 2.4.2 Optic Model Discretization

The discretization of the optic problem (16) is standard (see, e.g., [21, 25]). We consider the discrete space $W_h$ made of piece-wise affine continuous finite elements over $\Omega$, viz.,

$$W_h \stackrel{\text{def}}{=} \{w_h \in C^0(\bar{\Omega}) : w_h|_{K} \in \mathbb{P}_1 \quad \forall K \in \mathcal{T}_h\},$$

so that the fully discrete optical system reads as follows: For $1 \le n \le N$, find $\Phi_h^n \in W_h$ such that

$$\int_{\Omega} D^n \nabla \Phi_h^n \cdot \nabla w_h + \int_{\Omega} \mu_a^n \Phi_h^n w_h + \frac{1}{2} \int_{\Gamma_E} \Phi_h^n w_h + \frac{1}{6} \int_{\Gamma_S} \Phi_h^n w_h \\ + \frac{1}{2} \int_{\Gamma_T^e \cup \Gamma_B^e \cup \Gamma_R \cup \Gamma_D} \Phi_h^n w_h = 2 \int_{\Gamma_E} l(t_n) w_h \quad (20)$$

for all $w_h \in W_h$. Note that, at each time level $n$, the absorption and diffusion coefficients $\mu_a^n$ and $D^n$, depend on the results of the discrete poroelastic problem (19), namely on div $d_h^n$, via the relations

$$\mu_a^n \stackrel{\text{def}}{=} \left\{ \begin{array}{ll} (S\phi_0 + \text{div } d_h^n)\mu_{a,H_b O_2} + (1-S)\phi_0\mu_{a,H_b} & \\ + (1 - (\phi_0 + \text{div } d_h^n))\mu_{a,T} & \text{in } \Omega_T, \\ & \\ \mu_{a,B} & \text{in } \Omega_B, \end{array} \right. \quad (21) \\ D^n \stackrel{\text{def}}{=} \frac{1}{3(\mu_a^n + \mu_s^n)}$$

and with $\mu_s^n$ given by (15).
### 2.4.3 Discrete PPG Signal

Finally, the discrete radiant excitation flux trace on the receptor boundary $Γ_R$ can be evaluated via the following discrete counterpart of (17):

$$Q_R^n \stackrel{\text{def}}{=} \int_{\Gamma_R} \left( \frac{\Phi_h^n}{4} - \frac{D^n}{2} \frac{\partial \Phi_h^n}{\partial \mathbf{n}} \right). \qquad (22)$$

The complete solution procedure for PPG simulation is summarized in Algorithm 1. As discussed previously, there are two major steps in the PPG solver proposed in this paper. First, solve the poroelastic problem (19), to determine the local tissue deformation. Second, compute the fluence rate by solving the optical problem (20), with the tissue absorption coefficient accounting for the fluid volume fraction local variation via (21). Note that because the optical problem is quasi-static, one can solve it for a subset of all poromechanical time steps, to save computational time.

for $n = 1, \dots, N$ do

1. Solve poroelastic problem (19);

2. Update optical coefficients via (21);

3. Solve optical problem (20);

4. Compute total outflux on receptor via (22).

end

**Algorithm 1:** Numerical method for PPG signal simulation.

## 2.5 Experimental Framework

An experimental campaign has been devised to acquire PPG signals on the finger phalanx of a small set of healthy volunteers, with the purpose of validating and investigating the capabilities of Algorithm 1 for PPG simulation (see Section 3 below). Such a choice is to be understood as an initial step towards more challenging PPG body measurements locations. The corresponding experimental PPG dataset is publicly available at [16].
### 2.5.1 Acquisition Protocol

Measurements have been collected using the PPG sensor of the Withings BeamO (an OSRAM SFH7050 PPG module, see [39]), during an internal campaign at Withings with 20 healthy volunteers. The data collection process was performed as follows:

1. After informed consent, the subject's information is gathered;

2. The subject have to rest for 5 minutes;

3. A first arterial blood pressure (ABP) measurement is done in the left arm using Withings' blood pressure monitor (BPM Connect), an oscillometric sphygmo-manometer;

4. Subjects are told to rest their right index finger on the sensor and stay still. For each subject, 10 seconds of PPG is acquired twice;

5. Another ABP measurement is done in the same arm at the end to ensure the blood pressure stability during the procedure.

### 2.5.2 PPG Signal and Associated Biomarker

The PPG sensor used is composed of a set of three LEDs and a photodiode in a reflectance mode (see [39]). The LEDs are emitting at 3 different wavelengths in green, red and infrared domains (530, 655 and 940 nm, respectively). For each acquisition, we have a set of 3 channels of about 10 seconds of PPG signals, as depicted in Figure 3(a).

The PPG signal is frequently filtered in two components, the so-called AC and DC, with DC being the low frequency part and AC the pulsatile one, see Figure 3. More precisely, the DC component is extracted through a low-pass filter with a cut-off frequency, typically, between 0.1-0.5 Hz and AC through a pass-band filter with cut-off frequencies, frequently, between 0.5 and 10 Hz (see, e.g., [41] for a review on raw PPG filtering). Figure 3(b) shows the AC/DC ratio (the so-called perfusion index) with the associated biomarker considered in this work: the AC/DC peak-to-peak amplitude, per heart beat. There is a plethora of biomarkers in the literature to analyze the PPG
<figure><img src="image_5.png" /><figcaption>**Fig. 3** (a) Sample of the PPG signals acquired through the protocol described above, at the three wavelengths of interest, namely green, red and infrared; (b) Filtered PPG raw signals (from (a)) in its AC and DC components; (c) AC/DC signals, as presented in (b) with peak-to-peak AC/DC amplitudes, per heart beat.</figcaption></figure>
signal (see, e.g. [20] and the references therein). However, for the purpose of this work,
we focus solely on the AC/DC peak-to-peak amplitude, as no derivatives are required,
and it offers an initial assessment of the proposed modeling approach.

3 Results

The aim of this section is to investigate the performance of Algorithm 1 in the sim-
ulation of PPG signals. We first qualitatively asses Algorithm 1, in Section 3.1, by
illustrating how PPG measurements and simulations compare on a single heart beat
after calibration. In Section 3.2, we consider a broader set of model parameters with
the purpose of evaluating the model’s ability to reproduce the PPG data from the 20
measurements described in Section 2.5, by focusing on a specific biomarker (AC/DC
peak-to-peak amplitude). This analysis is a crucial step toward addressing the inverse
problem of estimating pulse pressure from PPG signals. In this context, we also exam-
ine whether, among the set parameter of parameters that yield simulation outputs
close to the PPG measurements, there exists at least one simulation whose pulse
pressure value approximates the measured brachial pulse pressure (mapped to the
finger).

All the simulations below have been carried out with an implementation of
Algorithm 1 in the FELiScE<fnref n="1" /> finite-element library.

**3.1 Qualitative Assessment of the PPG Simulation**

We start by presenting qualitative results of the behavior of Algorithm 1 on a single
example, for which we have measures of the PPG signal, for all three colors (green,
red, infrared), and of the brachial pressure, both systolic and diastolic (and thus of
the pulse pressure) as presented in Section 2.5.

<footnote n="1">¹https://gitlab.inria.fr/felisce/felisce</footnote>
The finger and PPG sensor geometries, together with the computational mesh, are presented in Figure 4. The geometry of the skin tissue (epidermis, dermis and hypodermis) comes from the Zygote<fnref n="2" /> library. The bone geometry has been segmented from the female cadaver of the Visible Human Project (see [1]). The sensor geometry was obtained from technical documentations of OSRAM's SFH7050 (see [39]). The LEDs and photodiode have a square geometry with sides lengths of 0.33 mm and 1.38 mm, respectively (see Figure 4(b)). All the geometries and meshes have been processed with the 3-matic<fnref n="3" /> software. The considered finger mesh is made of 163575 tetrahedra in the tissue region Ω<sub>T</sub>, and of 56385 tetrahedra in the bone region Ω<sub>B</sub>.

Poroelastic Simulation

In practice, since we employ a linear poroelastic model, the finger perfusion pressure P enforced in (19) corresponds to the difference between the instantaneous and diastolic pressure. This pressure waveform is derived from the ADAN model (see [54, 55]), from a left finger index distal artery, via the relation:

$$
P(t) = P_p \bar{P}(t), \qquad (23)
$$

where $P_p$ denotes pulse pressure (i.e., the difference of systolic and diastolic pressures) and $\bar{P}: \mathbb{R}^+ \rightarrow \mathbb{R}$ is the reference pressure profile shown in Figure 4(c).

The values of the poroelastic parameters are provided in Table 3. The simulations are performed for two cardiac cycles, with a time-step length of $\tau = 10^{-3}$ s. All the results presented below corresponds to the second cardiac cycle.

Figure 5(a) shows two clipped snapshots (lateral and frontal - across the PPG appa-
ratus) of the blood velocity magnitude in the deformed configuration (displacement
magnified by a factor 2), obtained with Algorithm 1, at the time instants $t = 0.08$ s

<footnote n="2">²https://www.zygote.com</footnote>

<footnote n="3">³https://wwwastic.com</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4 (a) Computational mesh with both a finger phalanx and the PPG sensor placement. Surface colors show boundary conditions labels; (b) OSRAM SFH7050 sensor geometry (see [39]); (c) Shifted and scaled left finger index distal artery pressure waveform $\bar{P}: \mathbb{R}^+ \to \mathbb{R}$, derived from the ADAN model (see [54, 55]).</figcaption></figure>

<table><caption>**Table 3** Summary of the parameters for the poroelastic model.</caption><thead><tr><th>Parameter</th><th>φ<sub>0</sub></th><th>k</th><th>ρ<sub>s</sub></th><th>λ<sub>s</sub></th><th>μ<sub>s</sub></th><th>ρ<sub>f</sub></th><th>λ<sub>f</sub></th><th>μ<sub>f</sub></th><th>P<sub>p</sub></th></tr><tr><th>Unit</th><th>%</th><th>cm<sup>2</sup> Ba<sup>-1</sup> s<sup>-1</sup></th><th>g cm<sup>-3</sup></th><th>Ba</th><th>Ba</th><th>g cm<sup>-3</sup></th><th>Ba s</th><th>Ba s</th><th>mmHg</th></tr></thead><tbody><tr><th>Value</th><td>0.082</td><td>0.012</td><td>1</td><td>1197213</td><td>49884</td><td>1</td><td>0</td><td>0.034</td><td>57</td></tr></tbody></table>

and $t = 0.32$ s, that is at maximal inflow and outflow, respectively; see also Figure 5(c). We clearly observe how the local volume of the finger increases as blood flows in. This is also consistent with the results of Figure 5(b), which provides the corresponding
<figure><img src="image_7.png" /><figcaption>**Fig. 5** (a) Clipped snapshots of the blood velocity magnitude in the deformed configuration (displacement magnified 2 times) at $t = 0.08$ s (maximal blood inflow velocity: top) and at $t = 0.32$ s (maximal blood outflow velocity: bottom), see Figure 4(c); (b) Clipped snapshots of the divergence of the displacement field at $t = 0.08$ s (top) and at $t = 0.32$ s (bottom); (c) Time histories, during the second cardiac cycle, of the (signed) maximal blood flow normal velocity on $\Gamma_T^e$, flow rate $Q$ and applied inlet blood pressure profile $P$.</figcaption></figure>

snapshots of the displacement divergence field, $\text{div } d$. Note that a positive divergence indicates an increase in the blood volume content. Finally, Figure 5(c) shows the time-histories of the (signed) maximal normal velocity on $\Gamma_T^e$, flow rate $Q \stackrel{\text{def}}{=} \int_{\Gamma_T^e} \mathbf{u} \cdot \mathbf{n}$ and pressure waveform $P$, given by (23). The maximal simulated inflow velocity attained here, slightly before the systolic peak, is about 23 cm/s. Such a value is in agreement with both measured values (see e.g., [7, Figure 3]) and simulated ones (see [54, 55] with a maximal value of about 16 cm/s in the left-hand finger index, assuming a parabolic profile).
Optical Simulation

The PPG geometry of the sensor platform is presented in Figure 4(b). The emitter boundary Γ<sub>E</sub> is made of the three components, Γ<sub>E</sub> = Γ<sub>E</sub><sup>G</sup> ∪ Γ<sub>E</sub><sup>R</sup> ∪ Γ<sub>E</sub><sup>I</sup>, respectively corresponding to the leds emitting in the green (530 nm), red (655 nm) and infrared (940 nm) wavelengths. The irradiant flux *l* on the emitter is assumed to be constant (in mW cm<sup>-2</sup>) on each activated LED, and given by

$$
l = \begin{cases} 16987 & \text{on } \Gamma_E^G, \\ 63542 & \text{on } \Gamma_E^R, \\ 42209 & \text{on } \Gamma_E^I. \end{cases}
$$

Such irradiance values are derived from the technical documentation of the LEDs and,
in particular, of the operating current of the LEDs. Note that only one led is activated
at a time, so that on the inactivated led we have *l* = 0.

<table><caption>**Table 4** Calibrated optical parameters.</caption>
  <thead>
    <tr>
      <th>Domains</th>
      <th colspan="6">Skin</th>
      <th colspan="2">Bone</th>
    </tr>
    <tr>
      <th>Parameters</th>
      <th>&epsilon;<sub>a,H<sub>b</sub>O<sub>2</sub></sub></th>
      <th>&epsilon;<sub>a,H<sub>b</sub></sub></th>
      <th>&mu;<sub>a,T</sub></th>
      <th>&mu;'<sub>s,T</sub></th>
      <th>S</th>
      <th>H</th>
      <th>&mu;<sub>a,B</sub></th>
      <th>&mu;'<sub>s,B</sub></th>
    </tr>
    <tr>
      <th>Units</th>
      <th>L cm<sup>-1</sup> mol<sup>-1</sup></th>
      <th>L cm<sup>-1</sup> mol<sup>-1</sup></th>
      <th>cm<sup>-1</sup></th>
      <th>cm<sup>-1</sup></th>
      <th>%</th>
      <th>%</th>
      <th>cm<sup>-1</sup></th>
      <th>cm<sup>-1</sup></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Green</th>
      <td>290</td>
      <td>300</td>
      <td>0.04</td>
      <td>10.0</td>
      <td>0.8</td>
      <td>0.397</td>
      <td>1.130</td>
      <td>22.24</td>
    </tr>
    <tr>
      <th>Red</th>
      <td>83.5</td>
      <td>800</td>
      <td>0.02</td>
      <td>1.3</td>
      <td>0.8</td>
      <td>0.397</td>
      <td>0.5287</td>
      <td>20.46</td>
    </tr>
    <tr>
      <th>Infra-red</th>
      <td>180</td>
      <td>200</td>
      <td>0.1</td>
      <td>1.0</td>
      <td>0.8</td>
      <td>0.397</td>
      <td>0.2210</td>
      <td>16.40</td>
    </tr>
  </tbody>
</table>

The AC and DC components of the measured PPG signals are used to (manually)
calibrate the parameters of the optical model (see Table 2), that AC and DC signals
obtained with Algorithm 1 match (for at least two significant digits) the measured
ones. Note that regarding the calibration process, on top of finding the absorption and
reduced scattering parameters, LED irradiance values were fixed within a ±25% range
<figure><img src="image_8.png" /><figcaption>**Fig. 6** (a) Snapshots of the fluence rate (logarithmic scale) obtained with Algorithm 1 for each color at *t* = 0.08 s, slicing the finger across each respective LED's color; as shown in (b); (b) Slices placement used in (a) along the LEDs' color of the PPG sensor, see also Figure 4(b), with green for green, red for red and black for infrared; (c) Comparison of the measured and simulated AC peak-to-peak (p-p) amplitude, DC, and (AC p-p amplitude)/DC; with G for Green, R for Red and IR for Infrared. (d) Measured and simulated AC/DC signals (after calibration on the measured DC and AC amplitude; see (c)).</figcaption></figure>

from manufacturer documentation to account for the LEDs irradiance uncertainty inherent from manufacturing variability. The resulting calibrated optical parameters are reported in Table 4.

Figure 6(b) provides snapshots within the finger of the fluence rate obtained with Algorithm 1 at *t* = 0.08 s, for each color, on frontal slices intersecting the LEDs, as shown in Figure 6(a). We can observe that the green penetrates less than red and infrared, which is consistent with observations commonly found in the literature (see, e.g., [3, 24]). Similar patterns are obtained at *t* = 0.32 s, which are omitted here
for the sake of conciseness. Figure 6(c) provides the comparison between measured and simulated PPG signals used for the calibration of the optical parameters, in terms of AC peak-to-peak (p-p) amplitude, DC, and (AC p-p amplitude)/DC. For the considered set of PPG data, the comparisons indicate that there exists a set of parameters so that simulated AC/DC amplitudes match the measured ones (up to two significant digits). Finally, Figure 6(d) shows the comparison of the measured and simulated AC/DC signal over a cardiac cycle (from the second cardiac cycle of the poromechanical simulation). Overall, the results show a good qualitative agreement, though significant deviations are clearly visible during the diastolic phase.

## 3.2 Simulated Optical Envelopes

In this section, we explore the capabilities of Algorithm 1 in reproducing the AC/DC peak-to-peak amplitude of the whole set of PPG measurements described in Section 2.5. The dataset comprises 20 measured cases, with about 10 heartbeats per case, yielding a total of 199 biomarkers for each color. The envelope of measurements hence consists of 199 points in $\mathbb{R}^3$ (one coordinate per color). The first goal of this section is to investigate whether the envelope of simulated biomarkers, provided by Algorithm 1, overlaps the measured one.

To this end, we use Algorithm 1 as a parametric solver and explore the parameters space (hypercube in $\mathbb{R}^{18}$) with a quasi-Monte Carlo method. We recall that the poroelastic solver (19) is decoupled from the optical solver (21), so that we can avoid sampling the whole parameters space at once. We hence begin by sampling the 5-dimensional poroelastic parameter space, defined by the bounds in Table 5, using a uniform distribution based on Sobol's sequences (see [46]). For each of these poroelastic simulations, the optical solver is invoked, for every color, with the parameters obtained from a uniform distribution (also using Sobol's sequences) in the 13-dimensional hypercube delimited by the bounds provided in Table 6.
## Parameters Ranges

As presented in Table 1, the poroelastic model involves 8 physical parameters. For the sake of simplicity, we assume that $ρ_f = ρ_s = 1$ and $λ_f = 0$ are fixed as in Section 3.1 (see Table 3). We express the Lamé coefficients in terms of the Young modulus and Poisson ratio, $E$ and $ν$ respectively. We assume here that the Poisson ratio remains constant, with $ν = 0.48$ (see [36]). Owing to [22, Table C1] and since we are interested in a single layer model for the skin tissues, we focus on bulk values for the Young modulus, ranging from $10^4$ to $10^6$ Ba. We recall here that the effective Young modulus is for the solid skeleton, and not exactly the linear elasticity Young modulus of the solid matrix. Such a point further reinforces the importance of parameters calibration here. It then remains three poroelastic parameters to characterize: $φ_0$, $k$, and $μ_f$. Based on various studies on blood volume content in the skin (see [30, 37, 42, 52]), we set the baseline volume fraction $φ_0$ of blood in the skin bulk to range from 1% to 40%. Due to the limited data available on the hydraulic conductivity $k$, we adopted an exploratory approach in which several poroelastic simulations were conducted across a wide range of $k$ values, and we then selected the range that produced physiologically plausible blood flow results (see Section 3.1). The determined range for $k$ is $10^{-3}$ to $1.6 \times 10^{-2} \text{ cm}^2 \text{ Ba}^{-1} \text{ s}^{-1}$. Blood viscosity depends on the hematocrit: viscosity increases with an increase in red blood cells. Here, following [27], we employ the 3-element Kriger-Dougherty model ([33]), that is

$$ \mu_f = \mu_m \left( 1 - \frac{H}{H^*} \right)^{-\eta}, $$

where $μ_m = 0.0123 \text{ Ba s}$ represents the plasma viscosity, $H^* = 0.99$ is the critical hematocrit volume fraction and $η = 2$ (see [40]). The hematocrit volume fraction $H$ is chosen to be within 0.36 and 0.5. Finally, the pulse pressure $P_p$ is chosen within the range 20 to 80 mmHg. These min-max values are inferred from linear regression on the
22 datasets of brachial and digital pulse pressures reported in [28, Table 2], which gives the following affine relation:

$$P_{\text{p,finger}} = 1.13 \cdot P_{\text{p,arm}} + 1.58, \quad r^2 = 0.63. \qquad (24)$$

In spite of the brachial to digital pulse pressure mapping being person specific, we nonetheless use (24) to estimate the finger pulse pressures of our set of 20 volunteers from their brachial measurements, yielding a minimal and maximal finger pulse pressure of 32 and 66 mmHg, respectively.

<table><caption>**Table 5** Summary of the poroelastic parameter ranges, bounding an hypercube in $\mathbb{R}^5$.</caption><thead><tr><th>Parameter</th><th>φ<sub>0</sub></th><th>k</th><th>E</th><th>H</th><th>P<sub>p</sub></th></tr><tr><th>Unit</th><th>%</th><th>cm<sup>2</sup> Ba<sup>-1</sup> s<sup>-1</sup></th><th>Ba</th><th>%</th><th>mmHg</th></tr></thead><tbody><tr><td>Min</td><td>0.01</td><td>10<sup>-3</sup></td><td>10<sup>4</sup></td><td>0.36</td><td>20</td></tr><tr><td>Max</td><td>0.4</td><td>1.6 · 10<sup>-2</sup></td><td>10<sup>6</sup></td><td>0.5</td><td>80</td></tr></tbody></table>

The poroelastic parameters ranges are summarized in Table 5.

As regards the parameters of the optical model (see Table 2), we assume here that the absorption and reduced scattering coefficients in the bone region, $\mu_{a,B}$ and $\mu'_{s,B}$ respectively, are fixed (one per color) as in example of Section 3.1 (see Table 4). The range of the hematocrit volume fraction $H$ is chosen as for the porolelastic solver in Table 5. For the arterio-venous saturation we assume $0.8 \le S \le 0.9$, which corresponds to an arterial oxygen saturation at 100% and a venous oxygen saturation between 0.6 and 0.8 (see [47]). Finally, for the rest of optical parameters in the tissue region $\Omega_T$, that is, $\epsilon_{a,H_bO_2}$, $\epsilon_{a,H_b}$, $\mu_{a,T}$, $\mu'_{s,T}$, there exists an extensive collection of tables and empirical relations with light wavelengths (see, e.g., [31]). However, these values are
<table><caption>**Table 6** Summary of the optic parameter ranges for each wavelength of interest, bounding an hypercube in $\mathbb{R}^{13}$.</caption><thead><tr><th>Parameters</th><th>ε<sub>a,HbO2</sub></th><th>ε<sub>a,Hb</sub></th><th>μ<sub>a,T</sub></th><th>μ'<sub>s,T</sub></th><th>S</th></tr><tr><th>Units</th><th>L cm<sup>-1</sup> mol<sup>-1</sup></th><th>L cm<sup>-1</sup> mol<sup>-1</sup></th><th>cm<sup>-1</sup></th><th>cm<sup>-1</sup></th><th>%</th></tr></thead><tbody><tr><td colspan="6">Green</td></tr><tr><td>Min</td><td>120</td><td>160</td><td>0.01</td><td>4.0</td><td>0.8</td></tr><tr><td>Max</td><td>480</td><td>640</td><td>0.1</td><td>16.0</td><td>0.9</td></tr><tr><td colspan="6">Red</td></tr><tr><td>Min</td><td>10</td><td>300</td><td>0.005</td><td>1.3</td><td>0.8</td></tr><tr><td>Max</td><td>150</td><td>1300</td><td>0.1</td><td>2.3</td><td>0.9</td></tr><tr><td colspan="6">Infrared</td></tr><tr><td>Min</td><td>30</td><td>50</td><td>0.03</td><td>1</td><td>0.8</td></tr><tr><td>Max</td><td>300</td><td>350</td><td>0.3</td><td>1.3</td><td>0.9</td></tr></tbody></table>

not utilized in the present study. This choice is consistent with our modeling assumptions, specifically the use of a homogeneous single-domain representation for the finger tissue $\Omega_T$. Since we lack precise information on the spatial distribution and individual concentrations of chromophores for a given subject, we do not attempt to model them explicitly. Instead, the parameter ranges (one per color) are centered around the calibrated values obtained from the reference example presented in Section 3.1, for which a physiological poroelastic simulation has been identified (see Table 4).

The optical parameters ranges are summarized in Table 6. Note that the saturation value $S$ in Table 6 is identical for all the three colors, reducing the number of independent parameters for the optical simulations to 13.

## Optical Envelopes

For this specific experiment, we consider 2400 samples from the poroelastic parameter space and generate 30 optical samples per color for each poroelastic case. This results in a total of 216000 simulations across the three color channels.

Figure 7 displays, as a scatter plot matrix, the AC/DC peak-to-peak amplitudes of both the 216000 simulations performed with Algorithm 1 and all the beat-to-beat
values of the 20 measured cases (which, as mentioned before, corresponds to an envelope of 199 measurement points in $R^3$). On the one hand, the pair red-infrared shows an excellent coverage, with coverage for all measured cases and every heartbeat. On the other hand, the pairs involving the green wavelength, while showing a good coverage, has some outliers (3 cases wholly out, that is all cardiac instances of these cases are not covered).

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Scatter plot matrix for all 3 colors delineating the envelopes of the simulations, alongside the measures of the AC/DC peak-to-peak amplitude biomarker.</figcaption></figure>
Now, an important question is whether, in a neighborhood of the measures (which will be precisely defined below), there exists a simulation whose pulse-pressure parameter $P_p$ is close to the measured value. By *measured value*, we mean here the brachial measured pulse pressure mapped to the finger using the relation (24). Since some green values are not covered by the simulated envelopes (see Figure 7), in the following we only consider the red and infrared outputs.

The notion of measure neighborhood is defined by taking into account the spreading of data. For each measured case $m$, $1 \le m \le 20$, we let $\mathcal{Y}^m = (y_{ij}^m)$ be the matrix of peak-to-peak AC/DC amplitudes, with $i = 1, 2$ indexing for the colors (1 for red and 2 for infrared) and $j = 1, \dots, n_b^m$ indexing for the heartbeat. The heartbeat mean is defined as the vector $\boldsymbol{\mu}^m = (\mu_i^m)$ in $\mathbb{R}^2$ with components

$$ \mu_i^m \stackrel{\text{def}}{=} \frac{1}{n_b^m} \sum_{j=1}^{n_b^m} y_{ij}^m, \quad i \in \{1, 2\}. $$

We can hence define the $2 \times 2$ measures' covariance matrix as

$$ \Gamma^m \stackrel{\text{def}}{=} (\mathcal{Y}^m - \boldsymbol{\mu}^m \otimes \mathbf{1}_{n_b^m}) (\mathcal{Y}^m - \boldsymbol{\mu}^m \otimes \mathbf{1}_{n_b^m})^\mathrm{T}, $$

where $\mathbf{1}_{n_b^m}$ stands for the array of ones in $\mathbb{R}^{n_b^m}$. We denote respectively by $\lambda_i^m$, $v_i^m$, with $i = 1, 2$, the eigenvalues and eigenvectors of $\Gamma^m$. Finally, we let $\{\boldsymbol{x}_l\}_{l=1,\dots,72 \cdot 10^3} \subset \mathbb{R}^2$ be the set of simulation outputs obtained with Algorithm 1. More precisely, the first and second components of $\boldsymbol{x}_l$ represent the red and infrared peak-to-peak AC/DC amplitudes, respectively, obtained from the $l$-th simulation. The simulated biomarker $\boldsymbol{x}_l$ is said to be within a neighborhood of the measured case $m$ whenever

$$ |\boldsymbol{x}_l \cdot \boldsymbol{v}_i^m| < 3\sqrt{\lambda_i^m} \quad \forall i \in \{1, 2\}. $$
An example of a case with 11 heartbeats, that is with 11 biomarkers per color, and
its associated simulation neighborhood is presented in Figure 8(a).

<figure><img src="image_10.png" /><figcaption>**Fig. 8** (a) Measured AC/DC peak-to-peak amplitudes biomarker for a single case with 11 heartbeats alongside simulations within the neighborhood of the measures. The “big-picture” is provided in Figure 7; (b) Bland-Altman plot of the mapped-to-finger brachially measured pulse pressure to the simulation nearest pulse pressure in a neighborhood of the red-infrared optical AC/DC peak-to-peak amplitude biomarkers (approximately 10 heartbeats per case) for 20 cases.</figcaption></figure>

Once the simulations within the neighborhood of the measures (per case *m*) are identified, we seek for the simulation that has the nearest pulse pressure (i.e., the parameter *P*<sub>p</sub> in (23)) to the mapped-to-finger measured pulse pressure (via relation (24)). The results for all 20 cases are presented in Figure 8(b). For these 20 cases, the mean of the difference in measured and simulation *P*<sub>p</sub> is 0.015 mmHg, with a standard deviation of 0.13 mmHg.

# 4 Discussion

In this work, we presented a mathematical model of photoplethysmography (PPG) applied to the finger, by combining a poroelastic model of tissue mechanics with an optical model of light propagation. The absorption coefficient within the optical model is modulated by the mechanical state of the tissue. In terms of validation, our objective was to provide a preliminary assessment of the model by simulating a range of parameter values and comparing the results to experimental data from a cohort of
20 volunteers, for whom both PPG signals and brachial blood pressure measurements were available.

We first presented results showing that, with one set of parameters (partially calibrated), the simulation provided by Algorithm 1 was able to reproduce both physiological literature results (e.g., blood flow velocity) and PPG measures (specifically, AC/DC peak-to-peak amplitudes). Encouraged by these promising initial results, we extended our evaluation by testing the proposed PPG simulator with a broader set of parameters on the small cohort of volunteers. Although the parameter intervals were heuristically set around the initial calibration, the results for AC/DC amplitudes were very good at the red and infrared wavelengths, but less consistent with the green. The rationale for such a loss of performance with the green is unclear. Principally, unlike red and infrared wavelengths, the penetration depth of the green wavelength does not reach the deeper dermal arterioles (see, e.g., [14]) hindering, for the green wavelength, the hypothesis that PPG originates from a variation in blood volume content, as described by (21) in our model with a single homogenized perfused tissue. In this regard, our results are consistent with those reported in [14]. As discussed, the parameters ranges in determining the envelopes were obtained heuristically, and expanding these ranges for the green wavelength might improve its coverage. However, this approach risks merely compensating for model inaccuracies specific to that wavelength. Therefore, adopting a multi-layer model to better represent the tissue structure could be a more effective strategy for enhancing green wavelength coverage. The proposed model can be readily adapted to incorporate a multi-layer perfused tissue, at the cost of a much finer mesh and an increase in the number of parameters.

Regarding the pulse pressure, we have shown that for all the measured cases in the red and infrared wavelengths, there exists a set of simulations (within the vicinity of the measured biomarkers) where at least one simulation's pulse pressure amplitude at the finger closely matches the mapped brachial-to-finger measured values, with a
standard deviation of less than 0.5 mmHg. This indicates that, within our cohort, we successfully identified a set of mechanical and optical parameters that reproduce both the optical outputs and pulse pressure measurements with high accuracy.

Our results were presented in terms of AC/DC peak-to-peak amplitudes, using a modulated (in amplitude) pulse pressure waveform from the literature (see Section 3.1). Significant differences between simulated and measured waveforms were observed, as shown in Figure 6(d). Clearly, incorporating additional biomarkers, notably temporal ones, would necessitate improving the temporal parameterization of the pressure waveform, for example, by accounting for heart rate variability or decomposing the waveform into a few Fourier modes.

Additionally, we have employed a single geometry for all simulations, despite the significant variation in finger sizes among individuals. Estimating the impact of a change in the shape of the finger with respect to a specific set of parameters is important and should be envisioned in further works. Also, we have assumed that the sensor is attached to the skin (i.e, mechanical contact was not modeled) and that the system is at equilibrium, so no variation in the sensor/finger force is accounted for, but it can be done within the linear poroelastic model. Such a point is important for instance when evaluating intra variability in the AC/DC peak-to-peak amplitude biomarkers, possibly caused by small forces variation, where AC/DC amplitude is expected to increase with small contact force (see, e.g., [48] and references therein).

In terms of perspectives, we currently see two major points for further studies. The first one regards the inverse and uncertainty quantification problems. Indeed, there, we have only scratched the issues and a significant amount of work is required for the assessment of the model error and robust parameters' estimation, in particular. The second one regards the use of the linear poroelastic model which limits the characterization of the blood pressure to its pulse pressure. While estimating pulse pressure has clinical significance (see, e.g., [17]), assessing the absolute pressure, especially systolic
and diastolic values, is highly desirable. However, to do so, the finite strain version of
the linear poroelastic model used here is required (see [12]).

Despite the limitations of the proposed optico-poroelastic PPG model, we are very encouraged by our initial results for PPG modeling and that it can be of interest, in particular, for the problem of cuffless continuous blood pressure monitoring and other challenges, such as the ones discussed recently in [13].

Declarations

**Ethics approval and consent to participate.** The PPG and pulse pressure data used in this study were collected by Withings SAS during an internal experimental campaign involving human participants. All participants provided informed consent prior to participation, and the study was conducted in accordance with relevant ethical guidelines and regulations. The data were anonymized before being shared for research purposes. As this paper involves only secondary use of anonymized data, no additional ethics approval was required.

**Consent for publication.** Not applicable. This work does not include any identifiable personal data or images.

**Availability of data and material.** The experimental dataset used in this paper is publicly available at [16].

**Competing interests.** The authors declare that they have no competing interests.

**Funding.** This work was financially supported by the "France Relance" recovery plan, launched by the French Government on December 21, 2020.

**Authors' contributions.** Conceptualization: AL, FV, MV, AB, DL, MF; Formal analysis and investigation: AL, JD, FV, MV, MW, AB, DL, MF; Funding acquisition: AB, DL, MF; Methodology: AL, FV, MV, DL, MF; Resources: AL, JD, MW, AB, DL,
MF; Software: AL, MF; Supervision: AB, DL, MF; Writing - original draft preparation: AL, MF; Writing - review and editing: AL, JD, FV, MV, MW, AB, DL, MF.

**Acknowledgements.** The authors would like to thank Dominique Chapelle, Céline Grandmont, Philippe Moireau and Romain Kirszenblat, for the fruitful and inspiring discussions on PPG modeling. We also express our gratitude to Marguerite Champion, for her valuable assistance with the computer implementation of the poroelastic solver in the FELiScE finite element library.

## Appendix A Generalized Inf-Sup Condition

The purpose of this section is to provide a proof the following generalized discrete inf-sup condition:

$$ \sup_{(v_h^s, v_h^f) \in \mathbf{V}_h^s \times \mathbf{V}_h^f} \frac{\int_{\Omega_T} q_h \operatorname{div}((1-\phi_0)\mathbf{v}_h^s + \phi_0 \mathbf{v}_h^f)}{(\|\mathbf{v}_h^s\|_{H^1}^2 + \|\mathbf{v}_h^f\|_{H^1}^2)^{\frac{1}{2}}} + \alpha (\mu_f s_h(q_h, q_h))^{\frac{1}{2}} \geq \beta \|q_h\|_{0,\Omega_T} \quad (A1) $$

for all $q_h \in Q_h$, and where $\alpha, \beta > 0$ denote two constants independent of $h$ and of $\phi_0$. This relation is crucial to ensuring that the discrete poroelastic problem (19) is well-posed. Indeed, using standard arguments (see, e.g., [21]) and (A1), one can show that (19) admits a unique solution for $\gamma_n$ sufficiently large.

The proof of (A1) proceeds as follows. We first note that from the numerical analysis of symmetric pressure stabilization methods for the Stokes problem there
holds (see, e.g., [9, 23]):

$$ \sup_{v_h \in V_h} \frac{\int_{\Omega_T} q_h \operatorname{div} v_h}{\|v_h\|_{H^1}} + (\mu_f s_h(q_h, q_h))^{\frac{1}{2}} \ge \beta^* \|q_h\|_{L^2} \quad \forall q_h \in Q_h, $$

with $\beta^* > 0$ a constant independent of $h$ and

$$ V_h \stackrel{\text{def}}{=} \{v_h \in [C^0(\bar{\Omega}_T)]^3 : v_h|_{K} \in [\mathbb{P}_1]^3 \quad \forall K \in \mathcal{T}_h^{\Omega_T}, \quad v_h|_{\Gamma_D \cup \Gamma_E \cup \Gamma_R \cup \Gamma_B \cup \Gamma_T^c} = 0 \}. $$

We hence have

$$ \sup_{v_h \in V_h^s} \frac{\int_{\Omega_T} q_h \operatorname{div} \left( (1 - \phi_0) v_h + \phi_0 v_h \right)}{\|v_h\|_{H^1} + \|v_h\|_{H^1}} + \frac{(\mu_f s_h(q_h, q_h))^{\frac{1}{2}}}{2} \ge \frac{\beta^*}{2} \|q_h\|_{L^2} \quad \forall q_h \in Q_h. $$

Therefore, since $V_h \subset V_h^s$ and $V_h \subset V_h^f$, we can set $v_h^s = v_h$ and $v_h^f = v_h$ in the previous relation, so that

$$ \sup_{(v_h^s, v_h^f) \in V_h^s \times V_h^f} \frac{\int_{\Omega_T} q_h \operatorname{div} \left( (1 - \phi_0) v_h^s + \phi_0 v_h^f \right)}{\|v_h^s\|_{H^1} + \|v_h^f\|_{H^1}} + \frac{(\mu_f s_h(q_h, q_h))^{\frac{1}{2}}}{2} \ge \frac{\beta}{2} \|q_h\|_{L^2} \quad \forall q_h \in Q_h. $$

The relation (A1) then follows by norm equivalence in the denominator, which completes the proof.

## References

[1] Ackerman MJ (1998) The visible human project. Proceedings of the IEEE 86(3):504–511. https://doi.org/10.3233/SHTI210988

[2] Allen J (2007) Photoplethysmography and its application in clinical physiological measurement. Physiological measurement 28(3):R1. https://doi.org/10.1088/0967-3334/28/3/R01
[3] Anderson RR, Parrish JA (1981) The optics of human skin. Journal of investiga-
tive dermatology 77(1):13–19. https://doi.org/10.1111/1523-1747.ep12479191

[4] Barnafi N, Zunino P, Dedè L, et al (2021) Mathematical analysis and numerical approximation of a general linearized poro-hyperelastic model. Computers & Mathematics with Applications 91:202–228. https://doi.org/10.1016/j.camwa.2020.07.025

[5] Barré M, Grandmont C, Moireau P (2023) Analysis of a linearized poromechanics model for incompressible and nearly incompressible materials. Evolution Equations and Control Theory 12(3):846–906. https://doi.org/10.3934/eect.2022053

[6] Biot MA (1962) Mechanics of deformation and acoustic propagation in porous media. Journal of applied physics 33(4):1482–1498. https://doi.org/10.1063/1.1728759

[7] Blair WF, Morecraft RJ, Brown TD, et al (1991) Transcutaneous blood flow measurements in arteries of the human hand. The Journal of hand surgery 16(1):169–175. https://doi.org/10.1016/S0363-5023(10)80037-1

[8] Brezzi F, Pitkäranta J (1984) On the stabilization of finite element approximations of the stokes equations. In: Efficient Solutions of Elliptic Systems: Proceedings of a GAMM-Seminar Kiel, January 27 to 29, 1984, Springer, pp 11–19, https://doi.org/10.1007/978-3-663-14169-3_2

[9] Burman E, Fernández M (2008/09) Galerkin finite element methods with sym-metric pressure stabilization for the transient Stokes equations: stability and convergence analysis. SIAM J Numer Anal 47(1):409–439
[10] Burman E, Fernández MA (2009) Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. Computer Methods in Applied Mechanics and Engineering 198(5-8):766-784. https://doi.org/10.1016/j.cma.2008.10.012

[11] Burtschell B, Moireau P, Chapelle D (2019) Numerical analysis for an energy-stable total discretization of a poromechanics model with inf-sup stability. Acta Mathematicae Applicatae Sinica, English Series 35:28-53. https://doi.org/10.1007/s10255-019-0799-5

[12] Chapelle D, Moireau P (2014) General coupling of porous flows and hyperelastic formulations—from thermodynamics principles to energy balance and compatible time schemes. European Journal of Mechanics-B/Fluids 46:82-96. https://doi.org/10.1016/j.euromechflu.2014.02.009

[13] Charlton PH, Allen J, Bailón R, et al (2023) The 2023 wearable photoplethysmography roadmap. Physiological measurement 44(11):111001. https://doi.org/10.1088/1361-6579/acead2

[14] Chatterjee S, Budidha K, Kyriacou PA (2020) Investigating the origin of photoplethysmography using a multiwavelength monte carlo model. Physiological measurement 41(8):084001. https://doi.org/10.1088/1361-6579/aba008

[15] Coussy O (2011) Mechanics and physics of porous solids. John Wiley & Sons, https://doi.org/10.1002/9780470710388

[16] Daraize J, Willemet M, Bendjoudi A (2025) Mathematical modeling of photoplethysmography: model assessment and validation. https://doi.org/10.57745/PCVIL2, Recherche Data Gouv, V1
[17] Dart AM, Kingwell BA (2001) Pulse pressure—a review of mechanisms and clinical relevance. Journal of the American College of Cardiology 37(4):975–984.
https://doi.org/10.1016/s0735-1097(01)01108-1

[18] Durian D (1998) The diffusion coefficient depends on absorption. Optics letters 23(19):1502–1504. https://doi.org/10.1364/OL.23.001502

[19] El-Hajj C, Kyriacou PA (2020) A review of machine learning techniques in photoplethysmography for the non-invasive cuff-less measurement of blood pressure. Biomedical Signal Processing and Control 58:101870. https://doi.org/10.1016/j.bspc.2020.101870

[20] Elgendi M, Liang Y, Ward R (2018) Toward generating more diagnostic features from photoplethysmogram waveforms. Diseases 6(1):20. https://doi.org/10.3390/diseases6010020

[21] Ern A, Guermond JL (2004) Theory and practice of finite elements. Springer, https://doi.org/10.1007/978-1-4757-4355-5

[22] Feng X, Li GY, Ramier A, et al (2022) In vivo stiffness measurement of epidermis, dermis, and hypodermis using broadband rayleigh-wave optical coherence elastography. Acta Biomaterialia 146:295–305. https://doi.org/10.1016/j.actbio.2022.04.030

[23] Fernández M (2013) Incremental displacement-correction schemes for incompressible fluid-structure interaction. Numer Math 123(1):21–65

[24] Finlayson L, Barnard IR, McMillan L, et al (2022) Depth penetration of light into skin as a function of wavelength from 200 to 1000 nm. Photochemistry and Photobiology 98(4):974–981. https://doi.org/10.1111/php.13550
[25] Groenhuis R, Ferwerda HA, Ten Bosch J (1983) Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. Applied optics 22(16):2456-2462. https://doi.org/10.1364/ao.22.002456

[26] Hansbo P, Hermansson J, Svedberg T (2004) Nitsche's method combined with space-time finite elements for ale fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering 193(39):4195-4206

[27] Hund SJ, Kameneva MV, Antaki JF (2017) A quasi-mechanistic mathematical representation for blood viscosity. Fluids 2(1):10. https://doi.org/10.3390/fluids2010010

[28] Imholz B, Langewouters GJ, Van Montfrans GA, et al (1993) Feasibility of ambulatory, continuous 24-hour finger arterial pressure recording. Hypertension 21(1):65-73. https://doi.org/10.1161/01.hyp.21.1.65

[29] Ishimaru A (1978) Wave propagation and scattering in random media, vol 1. Academic Press, https://doi.org/10.1016/B978-0-12-374701-3.X5001-7

[30] Jacques SL (1998) https://omlc.org/news/jan98/skinoptics.html/ [Accessed: 10/11/2023]

[31] Jacques SL (2013) Optical properties of biological tissues: a review. Physics in Medicine & Biology 58(11):R37. https://doi.org/10.1088/0031-9155/58/11/R37

[32] Jacques SL, Pogue BW (2008) Tutorial on diffuse light transport. Journal of biomedical optics 13(4):041302-041302. https://doi.org/10.1117/1.2967535

[33] Krieger IM, Dougherty TJ (1959) A mechanism for non-newtonian flow in suspensions of rigid spheres. Trans Soc Rheol 3(1):137-152. https://doi.org/10.1122/1.548848
[34] Kyriacou PA (2022) Introduction to photoplethysmography. In: Photoplethysmography. Elsevier, p 1–16, https://doi.org/10.1016/B978-0-12-823374-0.00001-3

[35] Kyriacou PA, Chatterjee S (2022) The origin of photoplethysmography. In: Photoplethysmography. Elsevier, p 17–43, https://doi.org/10.1016/B978-0-12-823374-0.00004-9

[36] Li C, Guan G, Reif R, et al (2012) Determining elastic properties of skin by measuring surface waves from an impulse mechanical stimulus using phase-sensitive optical coherence tomography. Journal of The Royal Society Interface 9(70):831–841. https://doi.org/10.1098/rsif.2011.0583

[37] Meglinski IV, Matcher SJ (2002) Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions. Physiological measurement 23(4):741. https://doi.org/10.1088/0967-3334/23/4/312

[38] Mejia-Mejia E, Allen J, Budidha K, et al (2022) Photoplethysmography signal processing and synthesis. In: Photoplethysmography. Elsevier, p 69–146, https://doi.org/10.1016/B978-0-12-823374-0.00015-3

[39] OSRAM (2016) Technical characteristics of LED SFH 7050 v1.1. Tech. rep., OSRAM Opto Semiconductors

[40] Ovarlez G, Bertrand F, Rodts S (2006) Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging. Journal of rheology 50(3):259–292. https://doi.org/10.1122/1.2188528

[41] Park J, Seok HS, Kim SS, et al (2022) Photoplethysmogram analysis and applications: an integrative review. Frontiers in Physiology 12:808451. https://doi.org/
10.3389/fphys.2021.808451

[42] Pasyk KA, Thomas SV, Hassett CA, et al (1989) Regional differences in capillary density of the normal human dermis. Plastic and reconstructive surgery 83(6):939-945

[43] Sarma PR (1990) Red cell indices. Clinical Methods: The History, Physical, and Laboratory Examinations 3rd edition

[44] Schmitt JM (1991) Simple photon diffusion analysis of the effects of multiple scattering on pulse oximetry. IEEE Transactions on Biomedical Engineering 38(12):1194-1203. https://doi.org/10.1109/10.137285

[45] Schweiger M, Arridge S, Hiraoka M, et al (1995) The finite element method for the propagation of light in scattering media: boundary and source conditions. Medical physics 22(11):1779-1792. https://doi.org/10.1118/1.597634

[46] Sobol' IM (1967) On the distribution of points in a cube and the approximate evaluation of integrals. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki 7(4):784-802. https://doi.org/10.1016/0041-5553(67)90144-9

[47] Solares G, Gomez-Fleitas M (2009) Central venous oxygen saturation: is it time to find its clinical applications? Journal of cardiothoracic and vascular anesthesia 23(2):272-273. https://doi.org/10.1053/j.jvca.2008.05.008

[48] Teng X, Zhang YT (2004) The effect of contacting force on photoplethysmographic signals. Physiological measurement 25(5):1323. https://doi.org/10.1088/0967-3334/25/5/020

[49] Ugryumova N, Matcher SJ, Attenburrow DP (2004) Measurement of bone mineral density via light scattering. Physics in Medicine & Biology 49(3):469. https://doi.org/10.1093/physrep/klh037
[50] Van Beekvelt MC, Colier WN, Wevers RA, et al (2001) Performance of near-infrared spectroscopy in measuring local o2 consumption and blood flow in skeletal muscle. Journal of applied physiology 90(2):511-519. https://doi.org/10.1152/jappl.2001.90.2.511

[51] Van Gemert M, Jacques SL, Sterenborg H, et al (1989) Skin optics. IEEE Transactions on biomedical engineering 36(12):1146-1154. https://doi.org/10.1109/10.42108

[52] Vidotto E, Koch T, Köppl T, et al (2019) Hybrid models for simulating blood flow in microvascular networks. Multiscale Modeling & Simulation 17(3):1076-1102. https://doi.org/10.1137/18M1228712

[53] Wang LV, Wu Hi (2007) Biomedical optics: principles and imaging. John Wiley & Sons, https://doi.org/10.1002/9780470177013

[54] Watanabe MSM (2013) Adan: um modelo anatomicamente detalhado da rede arterial humana para hemodinâmica computacional. PhD thesis, Laboratório Nacional de Computação Científica

[55] Watanabe SM, Blanco PJ, Feijóo RA (2013) Mathematical model of blood flow in an anatomically detailed arterial network of the arm. ESAIM: Mathematical Modelling and Numerical Analysis 47(4):961-985. https://doi.org/10.1051/m2an/2012053