Dynamic Modeling of Carbon Dioxide Transport through the Skin Using a Capnometry Wristband

The development of a capnometry wristband is of great interest for monitoring patients at home. We consider a new architecture in which a non-dispersive infrared (NDIR) optical measurement is located close to the skin surface and is combined with an open chamber principle with a continuous circulation of air flow in the collection cell. We propose a model for the temporal dynamics of the carbon dioxide exchange between the blood and the gas channel inside the device. The transport of carbon dioxide is modeled by convection–diffusion equations. We consider four compartments: blood, skin, the measurement cell and the collection cell. We introduce the state-space equations and the associated transition matrix associated with a Markovian model. We define an augmented system by combining a first-order autoregressive model describing the supply of carbon dioxide concentration in the blood compartment and its inertial resistance to change. We propose to use a Kalman filter to estimate the carbon dioxide concentration in the blood vessels recursively over time and thus monitor arterial carbon dioxide blood pressure in real time. Four performance factors with respect to the dynamic quantification of the CO2 blood concentration are considered, and a simulation is carried out based on data from a previous clinical study. These demonstrate the feasibility of such a technological concept.


Introduction
Respiration is a basic metabolic function of living organisms. On the one hand, it includes oxygen intake, which is used to produce energy at the cellular level. On the other hand, the CO 2 that is freed during the energy production is evacuated back into the environment. Most gas transportation within the organism is achieved via the circulatory system, whereas the gas exchange between the organism and the environment occurs at the level of the alveolocapillary barrier, although a secondary part of respiration takes place through the skin barrier [1].
The CO 2 in the blood is dissolved partly in its molecular species (aqueous carbon dioxide CO 2 (aq) and carbonic acid H 2 CO 3 ) and partly-due to the alkaline property of blood-in its deprotonated, ionic species (as carbonate CO 2− 3 or bicarbonate, which is also called hydrogen carbonate, HCO − 3 ). The proportion of CO 2 dissolved in plasma determines the pH of the blood [2]. The carbon dioxide pressure P CO 2 is linked to the concentration of the molecular species C CO 2 not = [CO 2 ] not = for "notation". Thus, an increase in the concentration of dissolved carbon dioxide induces a decrease in the blood pH. Metabolism is optimal for a pH of the blood in the range of 7.35 and 7.45, i.e., slightly alkaline [3]. It is the renal system, and the carbonic anhydrase enzyme in particular, that controls the blood's pH in the human body. An excess in pH might induce alkalosis and a deficit an acidosis. Both have a severe impact on human metabolism and the proper functioning of human organs [4]. Consequently, the monitoring of P CO 2 by capnometry is of major clinical interest.
Indeed, an increase in the carbon dioxide concentration might be a consequence of, for instance, the respiration track that is injured or obstructed, no longer evacuating CO 2 . This may not only happen in the case of chronic obstructive pulmonary disease (COPD) but also when an infectious disease, such as COVID-19, affects the lung. Capnometry is also useful to detect alveolar hypoventilation or hypercapnia. For clinical use, monitoring the carbon dioxide level is important for the follow-up of anesthesia or mechanically assisted respiration. It is used, among other things, for the monitoring of newborns in incubators. The main method to monitor respiration is to control the composition of the inhaled and exhaled air. An alternative method is to study gas exchange through the skin. When considering the development of a wearable device, this alternative solution is attractive. Controlling the air at the mouth or nose implies wearing a mask or a cannula, which is uncomfortable for longitudinal studies. Monitoring skin respiration-for instance, at the forearm-is more comfortable for the user. Developing a capnometry wristband would be of high benefit for homecare patient monitoring. However, it is a technological challenge to measure and monitor the low gas concentrations in a robust and reliable way. To this end, accurate instrumentation is necessary, in addition to relevant modeling and data processing. The model contributes to optimizing the design of the sensing device, by constructing the direct path on which the data processing will rely, and generating synthetic data to test the signal processing.
Several review papers on capnometry have been published [5][6][7][8][9][10]. The standard technology to measure carbon dioxide transcutaneous pressure is based on Severinghaus's electrochemical principle [11]. To obtain a shorter response time, we opt here for the optical non-dispersive infrared (NDIR) principle. This technology has been first introduced in the pioneering works of Thiele and Eletr et al. [12,13]. Some recent papers are proposing to use optical sensors based on fluorescent or luminescent dyes [10,14,15]. Currently, most of the capnometry devices rely on a model for a steady-state, static equilibrium between the carbon dioxide pressure within the blood and within the measurement device. In addition, most of the devices are based on a closed chamber measurement principle. In this work, we propose to consider an open chamber principle. Technology that combines, sequentially, an accumulation of carbon dioxide in a closed loop and a circulation of a nitrogen flow to flush the accumulated carbon dioxide has also been proposed [14,[16][17][18]. In the article of Iitani et al. [18], the accumulation phase lasted 60 s and the flushing lasted 30 s. The measurement is carried out at the beginning of the accumulation phase. The concentration is computed from the measurement of the accumulation flow rate. Each measurement is separated by the duration of the entire measurement cycle, which is one minute and thirty seconds.
In this work, we suggest using an open chamber principle including a continuous air flow. The use of nitrogen as carrier gas is not convenient for an autonomous wearable device. Hence, we propose to introduce air convection in the collection cell to speed up the carbon dioxide diffusion through the skin into the ambient air. However, this also comes with a dilution of the CO 2 concentration in the measurement cell. Hence, we need to study the compromise between the response time and the carbon dioxide gas concentration within the measurement cell. In order to study such an open chamber principle, we propose a model for the temporal dynamics of the carbon dioxide exchange between blood and the gas channel within the measurement device.
Signal processing methods include a large variety of approaches. An introduction is presented in the book chapter by Grangeat on signal processing from Nanoscience: Nanobiotechnology and Nanobiology [19], with a section dedicated to data extraction, including sub-sections on extracting physical quantities, systems approach, inverse problems, and regularized solutions, and sections on sensor correction and data analysis. In a pa-The IR radiation propagates along the measurement cell transversally with respect to the propagation direction of the carbon dioxide, which is normal to the skin surface, including multiple reflections on the lower grid in contact with the skin and the upper grid in contact with the collection cell. Two filters in front of each thermopile select the light wavelength to achieve the dual-wavelength measurement. We use a detection wavelength of 4.26 µm corresponding to a peak of the absorption intensity of carbon dioxide and a reference wavelength of 3.91 µm corresponding to a minimum of the absorption intensity for water. To obtain an autonomous wearable device, the electronics are embedded. They drive both the data acquisition, the pump, the heating, and the data exchange using Bluetooth Low Energy (BLE) protocols. Microelectromechanical system sensors (MEMS) are placed in each cell to measure continuously the temperature, pressure and relative humidity. The objective of this article is to describe a dynamic model of carbon dioxide transport through the skin on a capnometry wristband (named CAPNO) in order to introduce the theoretical framework for the development of this wearable device and to prove with simulated data the feasibility of this open chamber technology concept for the continuous monitoring of carbon dioxide blood concentration [43]. In Section 2, on Materials and Methods, we present the dynamic model that describes the CO 2 transport from bloodirrigated tissues through the skin to the measurement and the collection cells, followed by the methods for model evaluation. In Section 2.1, we present the continuous space and time dynamical model, including, in Section 2.1.1, a description of the four compartments' geometry (blood, skin, measurement and collection cells); in Section 2.1.2, a description of the measurement of the CO 2 concentration with dual-wavelength IR thermopiles; in Section 2.1.3, the assumptions leading to the boundary conditions; and, in Section 2.1.4, the physiological and physical parameters on which the model relies. Then, in Section 2.2, we describe the discrete state-space representation of CO 2 transport for the development of the simulation and inference software. This defines the direct model with the spatial discretization described in Section 2.2.1, the temporal discretization described in Section 2.2.2, the exogenous inputs described in Section 2.2.3, and the state-space model and the associated noise model described in Section 2.2.4. In Section 2.3, we present the model inversion to estimate the carbon dioxide concentration C in blood CO 2 and the associated pressure P in blood CO 2 at the input of the blood compartment from the carbon dioxide concentration C meas CO 2 computed from the NDIR optical measurement in the measurement cell. In Section 2.3.1, we derive a Kalman filter for the recursive estimation of the blood carbon dioxide concentration. In Section 2.3.2, we derive from the carbon dioxide blood concentration the carbon dioxide blood pressure. Then, in Section 2.4, we introduce the four performance parameters that we propose to characterize the quantification of the CO 2 blood concentration and its dynamic behavior. In Section 3, we present the results on simulated data, the performance parameters for the direct problem, the inverse problem and the complete model. We show also the performance parameters on a simulation of a realistic clinical test. In Section 4, we propose a discussion on these results and on the potential applications of this dynamic model. In Section 5, we present a general conclusion on the feasibility of this innovative capnometry wristband concept. We also discuss the contribution of the model to designing accurate instrumentation, to computing simulated data and to developing model-based signal processing.

Materials and Methods
The model is built upon the four compartments described in Figure 1. We consider a simplified one-dimensional (1D) diffusion-convection model of the transport equation along the z-axis perpendicular to the skin surface.

Continuous Space-Time Model of the CO 2 Transport
To study the transport of CO 2 species through the different media, we will start by writing down the conservation laws. Given a number N of moles of the species in a compartment of limited volume V, consider the net flux J (in mol·s −1 ) through its boundary Ω (according to the surface normals pointing outward from the volume), as illustrated in Figure 2. The transport equation is based on the conservation law of species in fluid mechanics [44]. The number of moles of the species in the volume changes with time according to the net flux J through the surface. Hence, writing down the number of moles within the volume V as a function of time leads to the following equation: The net flux through the surface Ω that governs the dynamics of the quantity N is an integration of local fluxes j (in mol·m −2 ·s −1 ) flowing through the boundary according to the surface normal → n (r) for all r ∈ Ω. We may thus write .
The transport equation is based on the conservation law of species in fluid mechanics [44]. The number of moles of the species in the volume changes with time according to the net flux through the surface. Hence, writing down the number of moles within the volume as a function of time leads to the following equation: where the minus sign stems from the fact that the surface normals are pointing outwards. If we consider an infinitesimal time step, we find that The net flux through the surface that governs the dynamics of the quantity is an integration of local fluxes (in mol·m −2 ·s −1 ) flowing through the boundary according to the surface normal ⃗( ) for all ∈ . We may thus write The fluxes ( , ) can be decomposed as a sum of three contributions: , which is given by Fick s first law of diffusion [45][46][47][48]: where the latter equality holds true only if the infinitesimal compartments are of constant volume (which is the case in Euclidean coordinates but is not the case in spherical nor in cylindrical coordinates). The diffusion constant ( ) depends on the medium and the species of .
, which is the media s convective flow, transporting the species with a speed ( ) in m·s −1 : Noting that a change in moles for a specific volume can now be expressed with the aid of the divergence theorem as This is a convection-diffusion equation that is the equivalent of Fick s second law [45,46,49] of a species concentration within a transporting medium that allows for a continuous concentration profile that is differentiable almost everywhere. In this case, we find the dynamics equations by considering that the left-and right-hand side must be valid for any subvolume; hence, the integrands must be equal in all coordinates ∈ : • Distributed source. If, in addition, we have a distributed source or sink within the volume that is, respectively, generating or absorbing moles at a constant rate, an additional term can be added to the right-hand side of convection-diffusion differential Equation (6) as a distributed source/sink term ( , ) (in mol⋅m −3 ⋅s −1 ). If the system interacts with the external world, modeled as exogeneous fluxes entering or leaving the system through the surface as a flux ( , ) (in mol⋅m −2 ⋅s −1 ), this can be added in the conservation equation as The fluxes j(r, t) can be decomposed as a sum of three contributions: , which is given by Fick's first law of diffusion [45][46][47][48]: where the latter equality holds true only if the infinitesimal compartments are of constant volume (which is the case in Euclidean coordinates but is not the case in spherical nor in cylindrical coordinates). The diffusion constant D(r) depends on the medium and the species of CO 2 .
• A convective flux j conv (r, t), which is the media's convective flow, transporting the species with a speed u(r) in m·s −1 : Noting that a change in moles for a specific volume can now be expressed with the aid of the divergence theorem as The transport equation is based on the conservation law of species in fluid mechanics [44]. The number of moles of the species in the volume changes with time according to the net flux through the surface. Hence, writing down the number of moles within the volume as a function of time leads to the following equation: where the minus sign stems from the fact that the surface normals are pointing outwards. If we consider an infinitesimal time step, we find that The net flux through the surface that governs the dynamics of the quantity is an integration of local fluxes (in mol·m −2 ·s −1 ) flowing through the boundary according to the surface normal ⃗( ) for all ∈ . We may thus write The fluxes ( , ) can be decomposed as a sum of three contributions: • A diffusive flux ( , ), which is given by Fick s first law of diffusion [45][46][47][48]: where the latter equality holds true only if the infinitesimal compartments are of constant volume (which is the case in Euclidean coordinates but is not the case in spherical nor in cylindrical coordinates). The diffusion constant ( ) depends on the medium and the species of .
which is the media s convective flow, transporting the species with a speed ( ) in m·s −1 : Noting that a change in moles for a specific volume can now be expressed with the aid of the divergence theorem as This is a convection-diffusion equation that is the equivalent of Fick s second law [45,46,49] of a species concentration within a transporting medium that allows for a continuous concentration profile that is differentiable almost everywhere. In this case, we find the dynamics equations by considering that the left-and right-hand side must be valid for any subvolume; hence, the integrands must be equal in all coordinates ∈ : • Distributed source. If, in addition, we have a distributed source or sink within the −j(r, t)· → n (r)dr = V ∇(D(r)∇C(r, t) − u(r)C(r, t))dr (6) This is a convection-diffusion equation that is the equivalent of Fick's second law [45,46,49] of a species' concentration within a transporting medium that allows for a continuous concentration profile that is differentiable almost everywhere. In this case, we find the dynamics' equations by considering that the left-and right-hand side must be valid for any subvolume; hence, the integrands must be equal in all coordinates r ∈ V: • Distributed source. If, in addition, we have a distributed source or sink within the volume that is, respectively, generating or absorbing moles at a constant rate, an additional term can be added to the right-hand side of convection-diffusion differential Equation (6) as a distributed source/sink term R(r, t) (in mol·m −3 ·s −1 ). If the system interacts with the external world, modeled as exogeneous fluxes entering or leaving the system through the surface as a flux j surf (r, t) (in mol·m −2 ·s −1 ), this can be added in the conservation equation The transport equation is based on the conservation law of species in fluid mechanics [44]. The number of moles of the species in the volume changes with time according to the net flux through the surface. Hence, writing down the number of moles within the volume as a function of time leads to the following equation: where the minus sign stems from the fact that the surface normals are pointing outwards. If we consider an infinitesimal time step, we find that The net flux through the surface that governs the dynamics of the quantity is an integration of local fluxes (in mol·m −2 ·s −1 ) flowing through the boundary according to the surface normal ⃗( ) for all ∈ . We may thus write The fluxes ( , ) can be decomposed as a sum of three contributions: • A diffusive flux ( , ), which is given by Fick s first law of diffusion [45][46][47][48]: where the latter equality holds true only if the infinitesimal compartments are of constant volume (which is the case in Euclidean coordinates but is not the case in spherical nor in −j(r, t)· → n (r) dr + Sensors 2023, 23,6096 The transport equation is based on the conservation law of species in fluid me [44]. The number of moles of the species in the volume changes with time accordin net flux through the surface. Hence, writing down the number of moles within ume as a function of time leads to the following equation: where the minus sign stems from the fact that the surface normals are pointing ou If we consider an infinitesimal time step, we find that The net flux through the surface that governs the dynamics of the quantity integration of local fluxes (in mol·m −2 ·s −1 ) flowing through the boundary accordin surface normal ⃗( ) for all ∈ . We may thus write where the latter equality holds true only if the infinitesimal compartments are of c volume (which is the case in Euclidean coordinates but is not the case in spherica −j surf (r, t)· → n (r) dr The concentration profile is not everywhere differentiable, this is especially so where the diffusion profile ( ) or the convection speed ( ) is discontinuous, which is the case at the interface between two different media. In this case, the continuity of the pressure and of the flux must be guaranteed at the interface, and the continuity equations yielding the appropriate interface boundary conditions are given as The concentration profile is not everywhere differentiable, this is especially so where the diffusion profile D(r) or the convection speed u(r) is discontinuous, which is the case Sensors 2023, 23, 6096 7 of 34 at the interface between two different media. In this case, the continuity of the pressure and of the flux must be guaranteed at the interface, and the continuity equations yielding the appropriate interface boundary conditions are given as where H − and H + are Henry's constants of the medium before and after the interface, respectively. The flux continuity implies that there is no build-up of specifies on the interface, the pressure continuity transcribes the fact that the partial pressure of CO 2 depends on its species in the specific media, and the pressure must be in equilibrium across the interface.
In what follows, the concentration of the specific species will label with a superscript referring to the medium or the compartment in which it is dissolved, considering the following four compartments: blood (blood-irrigated tissues), skin (dermal tissue), meas (measurement cell), and coll (collection cell). In addition, we will solely consider the transport equations according to an axis orthogonal to the skin surface (which will be labeled as the z-axis) considering that all net fluxes in the (x, y)-plane at any given zcoordinate are null (constant concentration in the plane, at least in the proximity of the considered axis), unless explicitly introduced.

Introducing the Four Compartments
When we refer to a compartment, we refer to a part of space in which the transporting medium is constant, which implies that the diffusion constant D is constant in that compartment. Four compartments will be distinguished: blood compartment, skin compartment, measurement cell and collection cell.
In the blood compartment, convection is determined by the net blood velocity, transporting the aqueous CO 2 .
The skin compartment defines a more rigid structure that does not allow for convection [45,46,50]. This is a purely diffusive medium.
The measurement cell has no forced convective air flow (the transporting medium) but houses the measurement device, which is a dual-wavelength IR thermopile, the functioning of which is detailed below.
The collection cell has a convective air flow in order to avoid any building up of CO 2 concentration in the measurement cell. However, higher convection speeds, although allowing for faster dynamics, will result in a reduction in the gas' concentration and hence a lower signal-to-noise ratio at the measurement device. In contrast, lower convection speeds will allow for build-up effects of the CO 2 concentration, resulting in an auto-regressive system favoring lower frequencies and, hence, slower dynamics.

Measuring CO 2 Concentration with Dual-Wavelength IR Thermopiles
The detectors of the measurement cell consist of two thermopiles sensitive to incident photon energy. A measurement and reference thermopile are used with their optical filters centered around the nominal wavelengths λ 1 = 4.26 µm and λ 2 = 3.91 µm, respectively. In order to reduce the contribution of other gases such as ambient air and water vapor, the NDIR differential measurement calculates the logarithm of the ratio of intensities at these two wavelengths.
U λ 1 and U λ 2 are the values read out from the thermopiles in the measurement cell in the presence of carbon dioxide, and U 0,λ 1 and U 0,λ 2 are the values read out from the thermopiles in the measurement cell in the absence of carbon dioxide. Considering the attenuation model of the Beer-Lambert's law, this attenuation depends linearly on the molar concentrations of the constituents. With the chosen wavelengths λ 1 and λ 2 , the attenuation measurement is a linear function of the molar concentration of carbon dioxide, independently of the concentration of water vapor. The attenuation profile k CO 2 of Sensors 2023, 23, 6096 8 of 34 CO 2 at the two wavelengths is very discriminating, i.e., k CO 2 (λ 1 ) k CO 2 (λ 2 ), but not discriminative for water ( In order to account for the non-linear effects associated with the wavelength spectrum of IR light, the multiple lengths of light pathways due to light reflection and the non-linear response of the thermopile IR sensors, in Grangeat et al. [51,52], a linear-quadratic model using a non-integer power of gas concentration has been described in order to relate the thermopiles' read-out voltage to the CO 2 concentration in the measurement cell. This model has been first proposed by Madrolle et al. for the quantification of a mixture of two diluted gases with a single metal oxide sensor (MOX) [53][54][55][56].
Knowing the logarithm of the voltage ratio, one can estimate the CO 2 molar concentration (C meas CO 2 ) in the measurement cell optical pathway using the following relationship (11) [53][54][55][56], where m, n, u are the (non-negative and real) quadratic model parameters.
If = ln U 0,λ 2 /U 0,λ 1 is not known, the intercept becomes another free parameter of the model: We place ourselves within the framework of a supervised calibration with the existence of N cal samples of known composition C meas CO 2 ,i . Using a Levenberg-Marquardt algorithm, the model parameters-m,n,u-are estimated by minimizing the following function Ψ:

Boundary Conditions for the System
The choice of boundary conditions is a prerequisite to solve the transport equations. In their article [57], Vaidya and Nitsche described the boundary conditions they have chosen for the simulation of the convection-diffusion equation of solutes in media with piecewise properties. In their article [36], Qazi et al. have described the boundary conditions they used to describe the kinetics of the carbon dioxide gas-liquid absorption on a membrane contactor.
At the interfaces, we will only conserve the diffusive flow, which means that the excess or default flow will leave or enter the system just before the interface.
For the direct model, we will consider the concentration at the extremal point of the known blood compartment (Dirichlet boundary condition). On the opposite side in the collection cell, we will suppose that there is no diffusive flow (zero derivative or von Neumann condition).
For the inverse model, we will consider the concentration obtained from the known thermopile voltage measurement (Dirichlet boundary condition) and suppose there is no diffusive flow at the extremal point of the blood compartment (von Neumann condition). Figure 2 illustrates the characteristic parameters used in the 1D transport model. Note that in order to allow a better visualization of the figure, the relative scales along the z-axis do not correspond to the relative size of each compartments.

Physical Parameters for the Transcutaneous CO 2 Transport
We have shown in this figure the inflow of blood. At the blood/skin interface, only the diffusive flux is conserved. We assume there is an outflow at the skin level within the blood to carry the convective flow, which does not cross the interface.
The air flow at the level of the collection cell allows mechanical convection to be established. We represent an incoming flow by the inflow arrows on both sides and an outgoing flow at the upper edge of the collection cell. The variables and parameters of the transport model of carbon dioxide from the blood to the ambient air through the skin and the measuring device are defined in the following tables. The parameters are taken from the literature.
Length and width parameters of Tables 1 and 2 are inherited from the device geometry. It defines the surface geometry, too.
Parameters defined in Tables 1-3 have been used to generate the data.

Physical Parameters for Blood Medium
The table's symbols are outlined in the following equations: where: • Q blood is the blood flow (cm 3 ·s −1 ); • A blood is the surface of the transverse section of the blood compartment (cm 2 ); • u blood z is the average blood velocity along the z axis (cm·s −1 ); • ∆x is the length of the blood compartment (cm); • ∆y is the width of the blood compartment (cm). Using the following formula, Formula (15), we compute the Henry constant for the blood based on the Ostwald's solubility coefficient for the blood: where: • R is the ideal gas constant-R = 0.0623637 m 3 ·mmHg·K −1 mol −1 ; • β blood is the Ostwald's solubility coefficient of the carbon dioxide in the blood at 42 • C (mol·m −3 ·mmHg −1 ); • T blood is the blood temperature (K).

Physical Parameters for the Skin
We consider hereafter that the skin compartment corresponds to the stratum corneum, which has the main effect on the diffusion of the carbon dioxide from the blood compartment to the measurement cell through the skin. Table 2 lists parameters used to characterize the skin. These are mean values, which are subject to intra-and inter-individual variations. For instance, skin elasticity might have an influence on skin compression induced by the wristband, which might have an influence on blood flow rate and average velocity.
Mass transfer coefficient k sc The diffusion coefficient D skin is computed from Equation (19): where: • α sc CO 2 is the Bunsen carbon dioxide solubility in the stratum corneum (mL(STPD)/mL solvent/mmHg) (STPD stands for standard temperature, pressure, and dryness); • K r sc CO 2 is the Krogh's diffusion constant (cm 2 ·s −1 ·mmHg −1 ). The Bunsen carbon dioxide solubility in the stratum corneum α sc CO 2 is given by the following expression: (20) where: • H skin is the Henry coefficient of the skin; • P air is the total air pressure at skin or membrane temperature.
If we suppose that H skin = 1.6 [33] and P air is 831.21 mmHg at 42 • C, we obtain: The Krogh's diffusion constant is given by the following expression: where: • k sc p is the mass transfer coefficient of the stratum corneum (m·s −1 ); • ∆z skin is the width of the stratum corneum (µm).
The mass transfer coefficient of the stratum corneum is computed from the conductance of the skin with respect to pressure difference per surface area using the following expression: where: • G vol di f f pres is the conductance of the skin with respect to pressure difference for the surface area of the device (mL·s −1 ·mmHg −1 ); • G vol di f f pres A blood is the volumic conductance of the skin with respect to pressure difference per surface area (mL·s −1 ·cm −2 ·mmHg −1 ); • A blood is the exchange surface area (cm 2 ).
The conductance of the skin is computed from the flow rate per exchange surface area coming out of the skin in the open air: where: • Φ out skin is the flow rate coming out of the skin (nl·min −1 ); • P blood CO 2 is the carbon dioxide pressure in the blood (mmHg); • P ambient air CO 2 is the carbon dioxide pressure in the ambient air (mmHg); • A blood is the exchange surface area (cm 2 ).
Itoh et al. [69] gives the following values for the skin conductance (volume flow rate) per surface area Typical value of CO 2 concentration in ambient air is 0.01613 µmol/cm 3 . We suppose here that a filter is placed in front of the air inlet of the collection cell, and this value is taken as null.

A Discrete Space-Time CO 2 Transportation Model
For an introduction to computer-aided modeling of material behavior, including physical and chemical parameters, and the mathematical tools for implementing those models in order to perform simulations, we refer to the referenced book of Jansens et al. [70] on computational materials engineering and, in particular, to chapter 5 by Kozeschnik on modeling solid-state diffusion [71]. The resolution of partial differential equations (PDE) using the finite difference method is described in the textbook by Langtangen and Ling [72]. Chapter 3 is dedicated to diffusion equations [73]. The textbook by Scherer on computational physics gives insight on numerical methods to compute physical models, in particular, for the discretization of differential equations [74].
We use a discrete state-space model in order to compute the time sequence of concentration profiles recursively along the transport path of the carbon dioxide, given a temporal sequence of carbon dioxide blood concentrations at the lower extremum of the blood compartment. This defines the recursive algorithm used for measurement data simulation; i.e., the previous time instant will provide the initial condition for the computation of the concentration profile along the transport path at the next time instant.
The following conventions for our notations are used: continuous scalar fields and constants will be denoted by light face characters. The space-time variables will appear as arguments between parentheses. Discrete representations will be given using arguments between square brackets; these quantities will appear in boldface lower-case characters when regrouped into a vector or boldface upper-case characters when regrouped into a matrix (a linear operator). Table 4 lists the parameters and notations used. In order to facilitate the discretization, we will use a finite difference scheme with regular, equidistant spatial sampling per compartment. At the boundaries we introduce a slack variable p interface as where the subscripts · − and · + refer, respectively, to the compartment below and above the interface (the z-axis pointing upwards). The slack variables can be eliminated from the equations using the interface conditions. As per Figure 3, we have N comp interior points for each compartment, together with two boundaries for which we have either a boundary or an interface condition, allowing to eliminate these from the set of unknowns. Given the thickness values of ∆z comp for the four compartments, we have the (relative) interface positions and may divide each compartment by placing N comp equidistant points, which are separated from their neighbours by a distance of δz comp = (N comp + 1) −1 ∆z. This results in a position vector, z = (z [1], z [2], · · · , z[n]) T , where the superscript · T is the transposition operator. If we note the concentration sampled in the ith point of the grid by c[i](t) = C(z i , t), we may form a concentration vector: We define a stencil for the first-order and second-order spatial derivative operators within a given compartment determining the ith equation using the following expression: The above system of equations then gives the following expression: We may thus express the spatial discretization of the convection-diffusion equation under matrix form: ∂ ∂t c(t) = F c(t) (29) where F is a linear finite-difference operator that encodes the spatial derivatives as well as the boundary and interface conditions defined below. Firstly, Dirichlet boundary condition is applied on the lower boundary. Then, the only equation that depends on the exogenous (boundary value) will be the one at index i = 1. We will have a new term D (δz) 2 + u 2δz c boundary (t) replacing the term that depends on c[0](t).
Secondly, Neumann boundary condition is imposed on the upper boundary (equation of index n). The condition states that there exists a point beyond the boundary such that And finally, Robin boundary condition is applied at each interface between model compartments. This boundary condition is imposed on a virtual interface point situated in between the ith and i + 1th sample points. We have the interface conditions that are built on one-sided differences in order to consider a single medium. At the interface, we have the limit concentrations in both media that are, respectively, given by c interface − and c interface + before and after the interface: From this, we achieve: This can be re-injected into the convection-diffusion equation, analogously: If the flow is purely diffusive through the interface, we have

Temporal Discretization
In order to integrate the above equation, we use an implicit Euler scheme [73] (chapter 3.2), [75,76] in which temporal sampling at regular intervals of δt leads to where Id n is the identity operator on R n . The discretized scheme is thus given by where the index k is the time sample index linked to the wall clock time, t 0 + (k − 1)δt and A def = (Id n − δt F) is the implicit linear operator. Since all eigenvalues of F are strictly negative, the matrix A has eigenvalues that are strictly bigger than one; hence, its inverse is a contractive operator, resulting in the stability of the algorithm. Note that to differentiate t (time) and k (time sample index), we use subscript for k and square brackets for [t].

Exogenous Inputs
To cope with exogeneous inputs in the system, such as the blood transport flow or the forced air flow through the collection cell, the equation is augmented with a steering matrix G and the exogenous inputs q that encode the boundary conditions as well as the flows: In discrete form, this gives G and q carry the information about the Dirichlet boundary condition at the lower bound, which gives and

State-Space Model and Noise
We define y k as the observation vector describing the measurement in the measurement cell at the time sample index k. As there is only one measurement, the vector y k is a scalar y k . The observation equation is where: • h is a canonical vector for which sole non-zero entry is associated with the position of the NDIR optical measurement (the last grid point in the measurement cell before the interface with the collection cell); • v k ∼ N 0, σ 2 v is the observation noise with a variance-covariance matrix R that describes the inaccuracies of sensor outputs as measurements are taken.
The difficulty in the inversion of the model is that no boundary condition can replace the knowledge of the concentration at the blood level. To overcome this difficulty, we suppose that the variations in concentration are slow and correspond to a first-order auto-regressive process, i.e., where: • w in blood [t] is the input noise model used to generate the input signal C(z blood,in , t); it is a white Gaussian noise process of variance σ 2 w . • ϕ is a signal regularity parameter of a first-order autoregressive signal model.
To derive the state-space equation, we combine this input signal model with the physical model describing the transport of the carbon dioxide. The associated state-space equations in discrete form read where δt·w k+1 ∼ N (0, Q) is the noise process that integrates the modeling errors. N (0, Q) is a normal vector distribution, with each component of zero average, and with a variancecovariance matrix Q. We assume that the modeling error noise w[t] and the observation noise v[t] are independent.

Kalman Filter Algorithm
The resolution of mathematical equations associated with diffusion phenomena is described in the reference book by Crank [77]. But in this book, the numerical methods do not consider noisy observations or state evolutions. Our problem is to estimate the state vector c k,k based on the previous state vector c k−1,k−1 and the (noisy) measurement y k . We propose to use a Kalman filter, which is an optimal filter that recursively computes a linear, least-mean square estimator [24][25][26]. The Kalman filter is adaptative in the sense that it tracks the noise level of both the measurements and the state evolution equations by updating the noise covariance matrices at each temporal step.
The Kalman algorithm works in a "prediction-correction" loop as described below and in [51,78]. For this description, we use the notations defined in Table 5.

Algorithm initialization:
During the initialization phase, at the time step k = 0, c 0,0 is initialized to null vector, and P 0,0 is initialized to identity matrix.

2.
Prediction step: Extrapolation (prediction) of the state and uncertainty covariance matrix at time step k + 1 from the current state (at time step k). We present below the equations for the implicit method used for time integration scheme. Equation (44) is for state extrapolation, and Equation (45)  We assume that the modeling error noise [ ] and the observation noise [ ] are independent.

Kalman Filter Algorithm
The resolution of mathematical equations associated with diffusion phenomena is described in the reference book by Crank [77]. But in this book, the numerical methods do not consider noisy observations or state evolutions. Our problem is to estimate the state vector , based on the previous state vector , and the (noisy) measurement . We propose to use a Kalman filter, which is an optimal filter that recursively computes a linear, least-mean square estimator [24][25][26]. The Kalman filter is adaptative in the sense that it tracks the noise level of both the measurements and the state evolution equations by updating the noise covariance matrices at each temporal step.
The Kalman algorithm works in a "prediction-correction" loop as described below and in [51,78]. For this description, we use the notations defined in Table 5. During the initialization phase, at the time step = 0, , is initialized to null vector, and , is initialized to identity matrix. 2. Prediction step: Extrapolation (prediction) of the state and uncertainty covariance matrix at time step + 1 from the current state (at time step ). We present below the equations for the implicit method used for time integration scheme. Equation (44) is for state extrapolation, and Equation (45) is for covariance extrapolation: 3. Correction step: State vector and uncertainty covariance matrix update using the estimates computed at time step − 1, , and , , and the current measurement . The Kalman gain given in Equation (46) expresses the weights given to the new noisy measurement or the estimated , state, which is also subject to different external disturbances: The estimation of the current state (a posteriori state estimator) of the system is performed by taking a linear combination between the estimate made at the previous A c k+1,k = c k,k + Gq k+1 ·δt (44) 3.
Correction step: State vector and uncertainty covariance matrix update using the estimates computed at time step k − 1, c k,k−1 and P k,k−1 , and the current measurement y k .
The Kalman gain given in Equation (46) expresses the weights given to the new noisy measurement y k or the estimated c k,k−1 state, which is also subject to different external disturbances: The estimation of the current state (a posteriori state estimator) of the system is performed by taking a linear combination between the estimate made at the previous moment k − 1 (a priori state estimator) and the new recorded data, as shown in Equation (47) below. Knowing the Kalman K k gain expression, we can give the expression to estimate the concentration c k,k at time k based on the difference between the current measurement y k and the estimated measurement hc k,k−1 from the previous state estimate. This difference defines the prediction error (y k − hc k,k−1 ), which evaluates the amount of new information brought by the current measurement: The estimation of the covariance matrix of the estimate error at time k from the estimated covariance matrix P k,k−1 at time k − 1 is based on Equation (48):

End of iteration:
The steps are iterated until either the end is decided by the user or the measurements are stopped. The output of the algorithm is the estimated CO 2 blood concentration corresponding to the first element of vector c k,k , at each time step.

Estimation of CO 2 Blood Pressure
Physicians are using CO 2 blood pressure P in blood CO 2 as a variable to characterize CO 2 blood concentration C in blood . This pressure is defined by the Henry law as the pressure the carbon dioxide would have in an air gas phase in equilibrium with the blood liquid phase: where β blood is the Ostwald solubility coefficient of the carbon dioxide in the blood. In this article, we are considering the carbon dioxide blood concentration variable.

Methods of Evaluation
The evaluation is conducted in three steps. Firstly, we evaluate the simulation of the direct transport problem. To model the CO 2 desorption for the four modelized compartments (blood, skin, measurement and collection), we used a minimum number of discretization points, 3 per compartment, plus the outer interface: N = 13. The following results are simulated using the parameters presented in Tables 1-3. The C in blood CO 2 input signal is simulated as a step transition between a normocapnia level and a hypercapnia level. The normocapnia level is defined as an average concentration (µ = 1.099 mol/m 3 ) of CO 2, which corresponds to a pressure of 40 mmHg.
Secondly, the results obtained for the transport inversion problem are illustrated on the previous simulations with different levels of noise (Gaussian zero mean white noise) added to the simulated observation vector to test the ability of the Kalman filter to adapt to the noise level.
Finally, we simulate a more realistic clinical test case simulating a sequence of stages with hypocapnia, normocapnia and hypercapnia levels. The Kalman inverse problem results are simulated.
The optimization of the numerical processing for the inversion problem is related to the quantification of the C in blood CO 2 CO 2 concentration level estimated in the blood medium. In these simulations, we assume that the CO 2 concentration in the incoming air is zero: C in col air CO 2 = 0. This is equivalent to assuming that a carbon dioxide filter has been placed over the air flow at the inlet of the collection cell.
The results are evaluated according to several performance factors with respect to the quantification of the CO 2 blood concentration estimated and to the time necessary for this quantification.
First performance parameter is calculated as the relative difference between the mean on the hypercapnia level (direct problem µ hypercapnia ) for input concentration C in blood CO 2 and the mean of the hypercapnia level (inverse problemμ hypercapnia ) of the estimated input concentration C in blood CO 2 . The hypercapnia means are calculated beyond the time when the level difference with the initial level (normocapnia) has reached 90% of the difference between the level at equilibrium after the transition (hypercapnia) and the initial level of normocapnia. In case of the "clinical" test, the evaluation is performed in the same way on the different successive capnia levels simulated (hypo, hyper1, hyper2 and normo capnia).
The second parameter is the root of the mean square error (RMSE) between the C in blood CO 2 concentration in the blood as input signal and the C in blood CO 2 estimated concentration in the blood as output signal. Computation of the RMSE is illustrated in Figure 4a. First signals are aligned, and then the root of the mean square error is calculated on the aligned signals. The rise time (t r ) corresponds to the difference between the time necessary for the initial level (normocapnia) to reach, respectively, 90% and 10% of the difference between the level at equilibrium after the transition (hypercapnia) and the initial level of normocapnia, and it is of interest for the temporal performance of the sensor. It is computed on the estimated concentration C in blood CO 2 . The method of rise time calculation is illustrated in Figure 4b. In the case of a "clinical" test, the rise time calculation is performed between the successive different capnia levels; indeed, some of them are fall time (normo to hypo or hyper to normo).
Other temporal performance parameters are the delays corresponding to direct, inverse and global problems. They are computed by analyzing the global temporal profile of the concentration along the transport of the carbon dioxide after placing a concentration step variation as input, as illustrated in Figure 4a.
The time difference corresponding to the maximum of the cross-correlation between the pulse at the input and the pulse at the level of the measurement cell defines the direct delay time (t d−dir ). This delay time is comparable to a propagation time of the carbon dioxide in the blood, the skin and the device up to the measurement cell. The time difference corresponding to the maximum of the cross-correlation between the pulse at the level of the measurement cell and the pulse at the level of the blood cell after the inverse problem defines the inverse delay time (t d−inv ). The time difference corresponding to the maximum of the cross-correlation between the pulse at the input and the pulse at the level of the blood cell after the inverse problem defines the global delay time (t d−global ).
These previous performance factors allow us to estimate two global performance factors. The relative performance (per f rel ) is defined as the ratio of the difference between the mean of the estimated and the input signal concentrations at the hypercapnia level.
The global signal-to-noise ratio rsb global is defined as the root mean square of the input signal and the RMSE: rsb global = 20× log10 rms C in blood where rms C in blood CO 2 is the quadratic norm of the input signal.

Evaluation of the Direct Approach
For all simulations, the input signal for our direct approach is illustrated in Figure 5. The input CO 2 concentration value C in blood CO 2 = 1.099 mol/m 3 is the average value (µ) of the concentration in the blood, which corresponds to a pressure of P in blood Hypercapnia level is defined as an increase of 10 mmHg, P hypercapnia CO 2 = 50 mmHg, which corresponds to a concentration of C hypercapnia CO 2 = 1.3738 mol/m 3 . Figure 5. System input signal. The variations in blood concentration are considered as input signal for the direct approach. There are two levels of blood concentration: one corresponding to a normocapnia level and one corresponding to a hypercapnia level.

Evaluation of the Complete Model for CO 2 Quantification
We consider three different noise variances added to the direct simulation result in order to see how they impact the time and to illustrate the adaptive property of the Kalman filter with respect to the noise level. We will also quantify the mean of the estimated hypercapnia blood level. We will examine the root mean square error between the real value of CO 2 blood concentration and the one estimated by the algorithm. The input signal of the Kalman filter C meas CO 2 is the result of the direct problem applied to a C in blood function with no noise, low noise and high noise added, presenting, respectively, a variance of 0 mol/m 3 2 , 10 −8 mol/m 3 2 and 10 −6 mol/m 3 2 , as illustrated in Figure 6.

Evaluation on a Realistic Simulated Clinical Case
In [41], Pierre Grangeat et al. described results of a clinical test conducted with a previous version of the CAPNO capnometry wristband in which the convection in the device channel was only induced by a difference in temperature between the ambient air and the heated skin, and the collection cell was under the measurement cell. We simulate such a test by considering a hypocapnia (decrease of 10 mmHg) phase, followed by a normocapnia phase and two different levels of hypercapnia (two successive increases of 5 mmHg), and then returning to the normocapnia condition. The objective of this test is to simulate a realistic experiment as it was conducted during a previous clinical evaluation. The simulation model is slightly different from the model described above. It is written using a compartmental model as described in [79].
The capnia levels and the duration of each phase are given in Table 6 and illustrated in Figure 7. We add two phases to the clinical chronogram:

•
An initialization phase to simulate the time that was used in the clinical test after installation of the device on the patient and before the test recording. During this initialization phase, we simulate a normocapnia level. • An extension phase in order to observe on simulations the return to normocapnia equilibrium. During this extension phase, we simulate a normocapnia level.  Figure 7 shows the capnia phase sequencing in the input of the model after an initialization phase in which equilibrium is reached and with an extension phase to observe the complete return to equilibrium. The hypocapnia phase lasts about 4 min, and the hypercapnia phases last about 6 min each.
We simulated realistic clinical data and examined the performance parameters.  The time necessary for the system initialization is visible in Figure 8a. For all other figures in the Results section, the initialization phase is omitted, as in Figure 8b. The figures are presented after the time necessary for the system initialization (i.e., after the time necessary to reach the normocapnia level as the system is initialized to 0).

Performance Parameters for the Direct Model
The optimization of the device time response involves the adjustment of certain parameters such as the ambient air flow entering the collection cell, which is responsible for the mechanical convection. This is necessary to extract the transcutaneous CO 2 flow, preventing the accumulation of CO 2 inside the device. An analysis related to the air flow traversing the collection cell is necessary. In a first intention, we will vary the ambient air flow between 10 −1 mL/min and 10 mL/min in order to obtain a partial pressure variation equal to the partial pressure associated with carbon dioxide in the ambient air.
The ambient air flow entering the collection cell contributes to the decrease in the response time of the system as well as a reduction in propagation time through the transport column but, in return, dilutes the concentration in the measurement and in the collection cell and, therefore, the level of the signal. There is a compromise to be made between lower time constants and the signal level at the output of the collection cell, as shown on Figure 9.

Performance Parameters for the Inverse Model
We tested the Kalman inversion algorithm on the result of the direct problem. We present the results successively for the three noise levels. In these results, the Kalman filter has been initialized to zero. The model noise variance is set to 10 −8 mol/m 3 2 , while the observation noise variance is set to 10 −6 mol/m 3 2 . The signal regularity parameter is ϕ = 0.

Noiseless Observation of a Step Function as Input
Firstly, we add no noise variance to the measured simulated output signal of the measurement cell. In Figure 10a, the Kalman filter estimation of the propagation of CO 2 through the different media in the four compartments is illustrated.

Low Noise Observation of a Step Function as Input
In the second test, we add a low noise variance to the measured simulated output signal of the measurement cell (cyan curve). Similarly, Figure 11a illustrates the propagation of CO 2 through the different media in the four compartments.

High Observation Noise of a Step Function as Input
Finally, in the third test, we add a high noise variance to the measured simulated output signal of the measurement cell (cyan curve). Similarly, Figure 12a illustrates the propagation of CO 2 through the different media in the four compartments.  Figure 12b. The estimated C in blood CO 2 curve is slightly noisy, but the mean hypercapnia level is well estimated.

Performance Parameters
The results obtained analyzing the direct and inverse approach are presented in Table 7. signal (curve in red) is illustrated in Figure 10b for the noiseless observation vector, in Figure 11b for low noise and in Figure 12b for high noise. We observe the delay to reach the hypercapnia measurement on each of these curves. For all noise levels tested, the hypercapnia level is reached.
These curves show the adaptability of the Kalman filter to react to the noise added to the signal. In all cases, the hypercapnia level is recovered. The Kalman filter has the advantage of acting as a real filter for the noise-to-signal ratio.
Delays and rise times are of the same order; they cannot be shorter than the physiology. Furthermore, the rise time of the global model is of the same order as the direct one, which makes our choice of a recursive Kalman filter as an inversion algorithm a realistic choice for a real-time device.

Performance Parameters for the Complete Model on a Realistic Clinical Test
Previously, we showed the adaptability of the Kalman filter for observation noise. Our purpose here is to show its adaptability to the model noise.

Compartment Model Simulation
The result of the simulation with the compartment model is illustrated in Figure 13. The durations of the different capnia levels are short; nevertheless, we observe the different levels on the observation curve (measurement cell). The Kalman inversion of the realistic clinical data simulation produces the result presented in Figure 14. This calculation is performed with (a) regularity parameter ϕ = 0 and (b) regularity parameter ϕ = −0.0036. Figure 15 presents detailed comparisons, firstly, in the blood cell, between blood CO 2 concentration in the input of the simulation model and the estimated CO 2 concentration using the Kalman filter (right scale), and, secondly (left scale), of the simulated measured CO 2 concentration in output of the simulation model with added noise (variance of 10 −6 mol/m 3 2 ) and the estimated CO 2 concentration in the measurement cell after Kalman inversion.  As we can see, the ϕ parameter allows us to compensate for concentration discrepancy in the blood cell. Magenta and red curves are on the same level in Figure 15b. However, we observe that the hypocapnia and hypercapnia levels are underestimated. The underestimation of the hypocapnia level is mainly due to the fact that the duration of the hypocapnia phase is too short.

Performance Parameters
The results are analyzed considering the six performance parameters defined above. They are presented in Table 8 depending on the parameters and presented in Table 9 for the global parameters, specifically. The main result of comparing the two columns of these tables according to the regularity parameter ϕ is that we see how the regularity parameter can act on the recovery of capnia levels. The relative performance parameter decreases from around 250% to less than 5%.
It shows also that, with this open chamber configuration (with forced convection), we can observe capnia phases of approximately 6 min time length, which is less than 10 min.
In [10], Dervieux et al. used a closed chamber model. They stated that if the sensor has a height of 1 mm, a 95% response will be achieved after 1 h 35 min and that, therefore, for a sensor to have a reasonable response time-e.g., below 10 min-it must be relatively thin, specifically in the range of 100 µm. We show in this work that, using an open chamber principle, this conclusion on the thickness of the device can still be released. Our simulated device has a thickness of 0.5 cm (measurement and collection cell thickness). This makes an open chamber device an alternative to a thin-film patch proposition.

Discussion
In this work, we propose an open chamber principle with a continuous circulation of air flow. But this requires the development of a dynamic model of carbon dioxide transport through the skin to recover the carbon dioxide blood content from the measurement sequence in the measurement cell. This allows the design of a model-based recursive signal processing approach based on a Kalman filter for a real-time estimation of the carbon dioxide blood pressure. A Kalman filter is relevant to estimate hidden variables in a noisy environment. The processing can be implemented with limited computational and memory resources.
This dynamic model described the transport of carbon dioxide from the blood to the collection cell through the skin and the measurement cell. It is based on convectiondiffusion equations in those compartments.
This model states all the parameters that have an influence on the measurement. These include parameters linked to the device architecture, to the IR source and sensors, to the operating mode, or to the patients. In Section 3, we illustrated the influence of the ambient air flow rate on the device time response and on the measurement level. Temperature also has an influence on both the physiological parameters of the model, such as blood velocity and blood flow, and on physical parameters such as the Henry coefficients, the solubility coefficients, the skin conductance and the associated mass transfer coefficient and Krogh's diffusion constant, the diffusion coefficients, and the carbon dioxide concentration and pressure in the ambient air. In this article, we have worked with a constant temperature of 42 • C. We have used mean values for skin parameters.
The variability of these parameters might have an influence on the measurement variability, including the intra-measurement variability linked to patient non-stationarity or measurement noise, inter-measurement variability for the same patient but for different operating modes or different devices, and inter-measurement variability among different patients.
The operating mode should also be well established. For instance, before starting the measurement, the device should be set on the patient, and some time should pass before starting the measurement in order to reach an initial equilibrium among the carbon dioxide contents in the four compartments.
The study of the measurement sensitivity to all these variabilities needs specific trials on experimental benches or on clinical trials. This should be performed to control the reliability of the measurement. For the parameters that are the more sensitive, complementary measurements should be achieved. This might be relevant, for instance, for patient-specific parameters including blood flow; skin thickness; skin vascularization; skin gas emission, such as water vapor induced by sweat or other volatile organic compounds (VOC); and contact between the device and the patient to control the skin heating or the gas loss according to the sealing of the carbon dioxide flow between the skin and the device. It might also be relevant for operating modes such as the variability of the ambient air parameters, including the carbon dioxide concentration, temperature pressure, humidity level, and flow speed. The measurement robustness is of primary interest for the development of a wearable device to be used in a homecare environment on moving people. The Kalman filter framework associated with the state-space model allows us to take into account the variability. In the simplest version, the variability is described by the two noise variables included in the model: the model noise, which might describe the errors linked by the use of an approximate state-space model, and the observation noise of the measurement error, including the inaccuracy error. The Kalman filter allows a recursive estimation of the variance of the noises and, thus, an improvement of the robustness with respect to measurement variability.
Parameters linked to the device architecture might be controlled by relevant calibration. The superiority of the IR measurement with respect to electro-chemical or electro-optical measurements is that frequent sensor calibration should not be required.
We have described here a simplified 1D version with a transport of the carbon dioxide gas along the axial direction perpendicular to the skin surface. The next step will be a two-dimensional (2D) version to take into account, using the same model, the transverse direction parallel to the skin surface and the axial direction. This 2D model will also allow us to consider a device architecture based on a transverse air convection flow parallel to the skin surface.

Conclusions
The development of a wearable device for the continuous monitoring of the blood carbon dioxide content is of primary interest for a homecare environment. It belongs to the large category of wearable health devices for vital sign monitoring described by Dias and Silva Cunha [80]. It belongs also to the general trend towards developing skin-based wearable devices as described by Jin et al. [81], wearable light sensors based on graphene material as described by Akinwande and Kireev [82,83], and wearable and miniaturized sensors for personalized and preventive medicine as described by Tricoli et al. [84]. The development of new room-temperature gas sensors based on a metal-organic framework nanocomposite as proposed by Zhang et al. [85] also gives a new perspective for using a gas sensor on a wristband. Tomasic et al. [86] have explained that the development of a wearable capnometry device is a requirement for the continuous monitoring of a COPD patient.
However, the development of such a capnometry device is a challenging topic since the transcutaneous carbon dioxide flow is very low: in the range of 0 to 800 nL/cm 2 /min. Therefore, we have proposed in this paper a dynamic model of carbon dioxide transport through the skin using a capnometry wristband. Such a dynamic model is a requirement not only for designing accurate instrumentation but also for computing simulated data and developing model-based signal processing. Thanks to this model, we have studied a new architecture in which the measurement cell is in contact with the skin and a continuous air flow is collecting the carbon dioxide gas. In this paper, we have proved with simulated data the feasibility of such a technological concept.
Such a model is a key contribution for the development of an autonomous wearable device as discussed in Section 4. The signal processing should be embedded to lower the data exchange rate between the wearable device and the smart phone or the computer that should store the data. It should also be resilient to measurement variability factors.
To improve the autonomy of the device, we have worked, in this paper, on the model to embed the data processing. But further studies will be required to lower the electrical power consumption. The main electrical power requirements are linked to data communication, skin heating and IR light emission.
We have described, in this paper, the transport of carbon dioxide from the blood to the collection cell. But such a model can be extended to other volatile molecular species, such as volatile organic compounds (VOC), including ethanol, acetone, isoprene or methane, and other diluted blood gases such as oxygen or hydrogen.
In [87], Mochalski et al. showed the possibility to measure other gases transcutaneously or to target other biomarkers. Other studies have been conducted, by Arakawa, to measure ethanol [88] or, by Ohkuwa, to establish a link between hypoxia and NO concentration [89]. In studying exposure to pollutants, Sekine et al. [90] established a relation between toluene emanating from the skin and inhalation exposure.