Using Drinking Data and Pharmacokinetic Modeling to Calibrate Transport Model and Blind Deconvolution Based Data Analysis Software for Transdermal Alcohol Biosensors (1)

Alcohol researchers/clinicians have two ways to collect subject/patient field data, standard-drink self-report and the breath analyzer, neither of which is passive or accurate because active subject participation is required. Transdermal alcohol sensors have been developed to measure transdermal alcohol concentration (TAC), but they are used primarily as abstinence monitors because converting TAC into more meaningful blood/breath alcohol concentration (BAC/BrAC) is difficult. In this paper, BAC/BrAC is estimated from TAC by first calibrating forward distributed parameter-based convolution models for ethanol transport from the blood through the skin using patient-collected drinking data for a single drinking episode and a nonlinear pharmacokinetic metabolic absorption/elimination model to estimate BAC. TAC and estimated BAC are then used to fit the forward convolution filter. Nonlinear least squares with adjoint-based gradient computation are used to fit both models. Calibration results are compared with those obtained using BAC/BrAC from alcohol challenges and from standard, linear, metabolic absorption, and zero order kinetics-based elimination models, by considering peak BAC, time of peak, and area under the BAC curve. Our models (with population parameters) could be included in a smart phone app that makes it convenient for the subject/patient to enter drinking data for a single episode in the field.


Introduction
Currently alcohol researchers and clinicians have only two ways to collect qualitative field data from research subjects and patients: the standard drink self-report and the breath analyzer. Neither of these data collection methods is passive in that they both require the active participation of the subject. This participation by the subject requires some level of training, is inconvenient, unnatural, and often results in missing, sparse, inaccurate, and imprecise data [42]. The primary pathway out of the body for alcohol consumed by humans is via an enzyme catalyzed metabolic chemical reaction in the liver [44]. However, smaller amounts of alcohol leave the body through other routes including respiration, urination and perspiration. Approximately 1% of the alcohol consumed by humans is excreted through the skin in the form of perspiration [39,45]. This transdermally excreted alcohol can be measured and has been shown to be correlated with BAC/BrAC [45]. Transdermal alcohol sensor (TAS) devices including the WrisTAS™ (Giner, Inc., Newton, MA) and SCRAM © (Alcohol Monitoring Systems, Inc., Denver, CO), have been developed to measure transdermal alcohol concentrations (TAC) [30]. These devices are currently used primarily as abstinence monitors (primarily for court ordered monitoring of individuals who have been convicted of driving under the influence) because of the difficulty of converting TAC data into more meaningful (to researchers and clinicians) BAC/BrAC data. This difficulty has two primary causes: 1) the complex and highly subject-dependent pathway that ethanol molecules follow from the blood through the skin to the TAC sensor, and 2) the current technological state-of-the-art with regard to the design and manufacture of a relatively lowcost and highly accurate electro-, or opto-chemical sensor. However, if, on the other hand, TAC sensor data could be processed and appropriately interpreted relative to BAC/BrAC, TAS devices have the potential to allow for the passive, accurate, and precise monitoring of BAC or BrAC over longer periods (i.e. for weeks at a time), and in fact, essentially continuously in time [11,45].
The requirement to have to first process TAC data before it can be used as a proxy for BAC/ BrAC is the result of the fact that these devices function like fuel cells by performing an oxidation reduction reaction that yields a count of ethanol molecules in the form of an electrical current. Consequently they measure alcohol rather than alcohol concentration. In addition, as we noted earlier, although they are bench calibrated, studies have shown that converting TAS measured TAC to estimates of BAC/BrAC can be highly device and/or subject dependent. This is in contrast to the breath analyzer which produces an estimate of BAC in the form of BrAC, by using an elementary stoichiometric result known as Henry's law [23] which describes the exchange of molecules between a liquid and gas that are in direct contact with one another. In this case, the underlying model involves a single parameter known as the partition coefficient which is the ratio of the concentrations of the molecule in the liquid and gas at equilibrium. It is generally accepted that this parameter varies only minimally from subject to subject and device to device. In the case of transdermal devices, on the other hand, when simultaneous BAC or BrAC and TAC are available, the TAC signal generally exhibits some degree of attenuation and latency when compared to the other two (we have used the qualifier generally here because on occasion the exact opposite has been observed!). Moreover, the level of attenuation and latency is often observed to vary from subject to subject, device to device and even drinking episode to drinking episode. Figure 1 above shows two plots of contemporaneous BrAC and TAC collected from two different subjects using two different devices. When the BrAC is compared to TAC, in the left panel, one can readily observe significant attenuation and latency, whereas in the right panel, there is virtually no attenuation at all and only minimal latency.
In earlier work [12,28,29,37,38], we have considered the problem of estimating BAC or BrAC from biosensor measurements of TAC based on the following general approach. We developed, analyzed, and made use of first principle mathematical models to create associated MATLAB software that converts TAC data into either BAC or BrAC estimates, depending on which kind of data was used to calibrate the models. More specifically, individualized physics-or physiologically-, based forward mathematical models in the form of partial differential equations with appropriate boundary conditions are fit to a particular subject and device using clinic or laboratory collected calibration data. The input to the model is BAC/BrAC and the output is TAC. Once the model is fit, it is used to invert the measured TAC data from drinking episodes in the field in the form of a deconvolution to obtain an estimate for the BAC/BrAC. The software uses two steps to obtain the BAC/BrAC estimates, which taken together can be viewed as a form of blind deconvolution. BAC/BrAC and TAC data from the laboratory calibration session are used to estimate a forward convolution filter and regularization parameters for a two or three parameter diffusion equation and its inversion. Then TAC data from field trial drinking episodes are used to deconvolve estimates, eBAC/eBrAC, of BAC/BrAC for the naturalistic drinking episodes in the field. Once an estimated BAC/BrAC curve has been obtained in this way, the software then produces estimates for three statistics that are of primary interest to alcohol researchers and clinicians: peak BAC/BrAC, time of peak BAC/BrAC, and area under the BAC/BrAC curve, along with the (essentially) continuous BAC /BrAC profile for the entire drinking episode. We first laid out this general approach in [12] and demonstrated the efficacy of our method on data collected in a clinical or laboratory setting. In [28] we reported on the development of a real-time assessment protocol that used a web application and a TAC sensor together with our BAC/BrAC estimator to assess changes in BrAC, alcohol responses, behaviors, and contexts over the course of a drinking event in thirty-two college students based on calibration data collected in a laboratory drinking session followed by a 2-week field trial. The functional analytic framework for the forward model we use in the present paper to describe the transport of ethanol from the blood through the skin to the TAC sensor in the form of a discrete-time parabolic PDE with unbounded (i.e. boundary) input and unbounded (i.e. point-wise) output along with the requisite approximation and convergence theory were developed in [29] and [37]. In addition, in those papers, we also developed and tested 1) the combined forward and inverse regularization procedure we employ here to guard against over-fitting when the forward model is calibrated and 2) a scheme to identify distinct drinking episodes in TAC data based on a Hodrick Prescott filter. Finally in [38] we looked at the modifications to our general theory that would be required so that various forms of boundary input (i.e. Dirichlet, Neumann, Robin, etc.) to the diffusion equation could be used in the forward model and how they affect the overall accuracy and performance of our data analysis scheme.
Unfortunately, at present, the need to collect calibration data in the clinic or laboratory before TAC data collected in the field can be processed and converted into equivalent BAC/ BrAC represents a significant obstacle to replacing either the self-report or breath analyzer with a TAC sensor. Indeed, supervising or participating in an alcohol challenge where the subject, wearing a TAC biosensor, is given an alcohol dose and blows into a breath analyzer at specified time intervals, represents a significant burden in the form of time and expense on the part of both the subject and the researcher or clinician. In this paper, we report on our efforts to greatly lessen this researcher, clinician and patient burden by eliminating the need to collect calibration data in the clinic or lab before the patient goes out in the field wearing only a TAC sensor. Our approach is based on replacing the alcohol challenge data with patient collected drinking diary data for the first drinking episode while the TAC sensor is worn. This can be done either in the clinic or lab (by giving the subject a measured dose of alcohol when they are fitted with the TAC sensor) or it can be done in the field by having the subject record their drinks for the first, or for that matter, any drinking episode. We then use this drinking data together with a nonlinear pharmacokinetic (Michaelis-Menten) based model for the absorption, metabolism, and elimination of alcohol from the bloodstream to obtain an estimate of BAC/BrAC for this drinking episode. This estimate of BAC/BrAC together with the TAC for this episode is then used to calibrate the forward model that describes the transdermal transport of ethanol from the blood to the TAS sensor and its processing and analysis once it gets there. This forward model is then used as a basis for the deconvolution of BAC/BrAC from the TAC data collected in the field by the TAS sensor for all subsequent drinking episodes that occur while the sensor is being worn. Note that no additional active subject input in the form of a drinking diary or any other voluntary compliance other than wearing the sensor is required for the estimation of the BAC/BrAC that results from these subsequent field drinking episodes. We note further that the drinking diary for the first drinking episode that we do require is of relatively low resolution. That is, we simply require the number of standard drinks consumed and the times at which they were served during the episode. A standard drink as defined by the National Institute of Alcohol Abuse and Alcoholism (NIAAA), one of the research institutes of the NIH, is one that contains 14 grams or 0.56 fl. oz. of pure ethanol [8,21]. For example, any of the following constitute one standard drink: one regular 12 fl. oz. beer, 5 fl. oz. of table wine, 2-3 fl. Oz. of aperitif, or 1.5 fl. Oz. or a shot of 80-proof spirits (hard liquor). In addition, our scheme is such that if a higher resolution drinking diary is available (as it is likely to be in the case of a controlled study involving paid subjects or volunteers) we can make use of this more complete dataset to obtain a more accurate pharmacokinetic estimate of BAC or BrAC. Our approach is illustrated schematically in Figure 2. Note that we label Stage 1 as being off-line in the sense that the parameters in the drinking/elimination model used to produce estimated BAC/BrAC would be fit based on stratified (with respect to body mass, sex, drinking history, etc.) population data. We also consider Stage 3 to be quasi-offline in the sense that it would only be used for the first drinking episode, the one for which estimated BAC/BrAC is available.
Other authors have also looked at the problem of extracting BAC/BrAC like data from TAC measurements. In 2007 and 2008 papers [51,52] Webster and Gabler looked at a coupled Michelis-Menten based elimination model and diffusion model for transport through the skin in the context of designing an auto ignition interlock system to prevent driving under the influence. However they only considered the forward simulation problem and were primarily interested in studying and characterizing the observed lag between TAC and BAC/ BrAC and in particular how it depended on covariates such as gender, body mass, and amount of alcohol ingested. They used values gleaned from the literature for the physical and biophysical parameters that appear in their model rather than trying to obtain them from data in the form of an inverse problem. In a sequence of papers in 2012, 2013, 2014 and 2015, [11,17,18,19] Dougherty and his group used standard linear regression to fit finite dimensional models for peak BAC and time to peak BAC and total number of drinks consumed to factors such as peak TAC, time to peak TAC, area under TAC curve, and sex, etc. and to then use the resulting fits to draw conclusions about latency, attenuation and dependence.
An outline of the remainder of the paper is as follows. In Section 2 we describe our transdermal alcohol transport model, the corresponding blind deconvolution problem, our two-stage process for first calibrating the model or equivalently identifying the convolution filter (Stage 3 in Figure 2 above), and then how we use the resulting calibrated model to deconvolve BAC/BrAC from field TAC measurements (Stage 4 in Figure 2 above). In Section 3 we consider our nonlinear pharmacokinetic based model that takes the number of standard drinks ingested as input and produces estimated BAC/BrAC as output (Stage 2 in Figure 2 above), and we then describe the process by which this model would be fit to stratified population drinking data/BrAC data (Stage 1 in Figure 2 above). In Section 4 we present our numerical results for fitting the nonlinear model to data. In Section 5 we compare the use of our nonlinear pharmacokinetic model for obtaining calibration data to the use of a number of existing linear/zero order kinetics absorption/elimination alcohol models and to the use of alcohol challenge data in calibrating our transdermal alcohol model as part of our TAC sensor data analysis software. Section 6 contains a discussion of our findings and a few concluding remarks.

BrAC from TAC
We seek a linear input/output model for TAC, y, as a function of BAC/BrAC, u, in the form of a convolution [9] (2.1) where the unknown convolution filter or kernel, k, is parameterized by a finite dimensional vector of parameters, q ∈ Q, Q the set of feasible parameters. Moreover, to keep the dimension of the parameter space as low as possible, rather than simply parameterize k in some arbitrary fashion (for example, as a linear combination of a set of basis functions such as orthogonal polynomials or splines), we take advantage of basic physical and physiological principles that govern the transport of ethanol from the blood through the skin. Thus we consider a first principles forward model for the transdermal transport of alcohol based on the assumption that the ethanol molecules diffuse through the interstitial fluid that surrounds the cells that make up the epidermal layer of the skin (see also [2,12,49,52]). The resulting system of equations involve parameters that essentially describe the rate at which ethanol diffuses through the epidermis, the effective net rate at which ethanol enters and leaves the skin, and how it is processed by the transdermal alcohol sensor, etc.
Letting φ(t,x) denote the concentration of ethanol in the epidermal interstitial fluid in units of mg/dl at depth x and time t, our basic underlying model in the form of an initial-boundary value problem for a one dimensional diffusion equation with a Robin boundary condition on one end and a Neumann boundary input on the other. We note that other modeling choices for the boundary conditions would also be reasonable [38]. The model takes the form (see [37] for the details) where the input, u, is BAC/BrAC, the output, y, is TAC, both in units of mg/dl, T, in units of hours (h), is arbitrary with σ > 0, and ρ,D,L, α,β, and γare the usual physical parameters associated with diffusion that we will use to fit the model to observations of u and y. It is not surprising that not all six parameters are independent. In fact (see [37]), by making the change to the dimensionless independent variables and , setting T = 1, choosing an appropriately weighted (the weight is a function of the parameters) L 2 inner product when we convert the system to weak form, and taking advantage of linearity and the zero initial conditions, we find that without any loss of generality we may consider the simplified transformed system where q 1 and q 2 are two unknown dimensionless and, in general, subject-and devicedependent parameters. Note that our basic modeling assumptions are that the ethanol flux through the upper surface of the skin (x = 0) is proportional to the concentration of ethanol at the surface and that the ethanol flux from the dermal layer of the skin into the lower surface of the epidermal layer (x = 1) is proportional to the concentration of alcohol in the blood (or breath as a surrogate). We have assumed that there is no alcohol in the epidermal layer at time t = 0 (we will deconvolve one drinking episode at a time with the beginning of the episode defined to be when the BAC/BrAC rises above zero and the end of the episode being when the BAC/BrAC and TAC fall below a pre-specified tolerance and remain there for a pre-specified period of time). We assume that the TAC sensor measures alcohol concentration (i.e. TAC) at the upper surface of the epidermal layer (x = 0) . We note that in this context, we have also used a three parameter hyperbolic diffusion model which exhibits finite speed of propagation [34] and obtained similar results to the ones we present here.
To obtain y in the form of u convolved with the convolution filter or kernel k(·) = k(q, ·) as in (2.1), we reformulate (2.2) -(2.5) in weak or variational form [27] and then rely on the theory of linear semigroups of operators [20,22,35]. Let H = L 2 (0,1) together with the standard inner product , and norm denoted by | · |, and let V be the Sobolev space [1] V = H 1 (0,1) with norm denoted by ∥ · ∥. Then we have the usual dense and continuous embeddings , where V* denotes the space of distributions dual to V [27,41,46,53]. For q = (q 1 , q 2 ) ∈ Q, a compact subset of R + × R + , we define the bilinear form a(q; · , ·) : V × V → R and the functions b(q; ·) : V → R and c(·) : V → R by where 〈 ·, · 〉 denotes the natural extension of the H inner product to the duality pairing between V and V*. Note that because our model (2.2)-(2.5) involves boundary input and pointwise observation, also on the boundary, special care will have to be taken in recasting the system in terms of an abstract evolution equation using semigroups. Indeed, the resulting input and output operators turn out to be unbounded with respect to the standard state space (L 2 (0,1)) where one would typically set the problem [5,10,15,36,43,47].
For q ∈ Q, the q-dependent bilinear form on V × V, a(q; ·, ·) : V × V → R defined in (2.6), defines a bounded linear operator by , for ψ 1 ,ψ 2 ∈ V. Then, letting denote any of the Hilbert spaces V, H, or V*, we consider the linear operator given by A(q)ψ = A(q)ψ for . It can be shown [5] that A(q) is a closed, densely defined unbounded linear operator and the infinitesimal generator of an analytic semigroup of bounded linear operators, on .
We write the model equations (2.2) -(2.5) or, equivalently, (2.9) -(2.11) in strong form as an abstract evolution equation initial value problem in V* as with output equation Then using the fact that {e A(q)t : t ≥ 0} is an analytic semigroup on V*and therefore that , for ψ ∈ V*, we obtain from the abstract variation of constants formula It follows that for q = (q 1 ,q 2 ) ∈ Q, y is given by (2.1) with the convolution filter k(·) = k(q; ·) given by a well-defined real valued function of t defined for all values of t > 0. Now let a sampling time or interval τ > 0 be given and set φ i = φ(iτ, ·), y i = y(iτ) and u i = u(iτ), i = 0,1,2,... Then, defining the zero-order hold input u(t) = u i , t ∈ [iτ,(i + 1)τ), i = 0,1,2, ... , the variation of parameters formula yields where and . Boundedness of the operators Â(q) and B(q) follows from the fact that {e A(q)t :t ≥ 0} is an analytic semigroup on V,H and V* [5]. Moreover, it can be argued (see [37] defined with respect to the usual uniform mesh, , and set (note the are the usual "pup tent" or "chapeau" functions of height one and support of width 2/n, given by ). Let P n :H → H n denote the orthogonal projection of H onto H n with respect to the H inner product. It is not difficult to argue that in H for ψ ∈ H and V for ψ ∈ V [3,4]. For n = 1,2,..., and q ∈ Q, define to be the finite dimensional linear operator , where is the orthogonal projectino of V onto H n with respect to the inner product 〈·,·〉 a = a(q; ·,·) on V. For i,j = 0,1,2, ..., n, is its matrix representation. Finally we set Â n (q) = e A n (q)τ and nothing that b(q) ∈ H n for all n = 1,2, ..., we obtain i = 1,2, ..., where b̂n(q) denotes the vector representation of b(q) ∈ H n with respect to the basis . The approximating discrete time convolution kernel is then given by Letting {ũ j ,ỹ j }, j = 0,1,2, ..., N denote calibration data in the form of either BAC-TAC or BrAC-TAC pairs for either a clinic or laboratory alcohol challenge session or the first drinking episode, we formulate the calibration problem as a constrained least squares fit to data. We seek in Q which minimizes: We denote the calibrated convolution kernel by The optimization of (2.13) is achieved using an iterative constrained gradient based search.
We compute the gradient using the adjoint method [24,37]; where for i = 0,1,2,..., N, satisfies the adjoint system with . We note that the gradient can be computed directly. The tensor can be computed at the same time that Â n (q) is computed using the sensitivity equations [32]. For t ≥ 0 and q ∈ Q. set Φ N (q;t) = e A N (q)t . Then Φ N (q;·) is the unique principal fundamental matrix solution to the initial value problem Then, differentiating with respect to q, interchanging the order of differentiation, using the product rule, and setting we find that

It then follows that
Once calibration is complete and the parameters and therefore the convolution kernel, , i = 1,2,...., as well, have been estimated we formulate the deconvolution problem for field BAC/BrAC as a linear least squares fit to data with an added positivity constraint [26]. (2.14) where , j = 0,1,2,...,ν. Our optimal estimate for the BAC/BrAC signal û(t) corresponding to TAC signal ŷ(t) is then given by .
Finally, we note that to mitigate the effects of over fitting due to the inherent ill-posedness of the deconvolution problem, terms quadratic in the U i ′s are weighted and added to the expression for J D given in equation (2.14). The performance index for the deconvolution now becomes (2.15) where J D is given by equation (2.14), R 1 and R 2 are (M+1)×(M+1) Gramian matrices for and , respectively (i.e. and ), and r 1 and r 2 are non-negative regularization weights. For given nonnegative values of the regularization weights, r 1 and r 2 , finding , j = 1,2,...,M that minimize to equation (2.15) is a straight forward, linear-quadratic programming problem for which accurate and efficient algorithms and software are readily available (e.g., the MATLAB routine LSQNONNEG) [6,7]. To estimate the regularization weights we further modify the calibration procedure and add a second calibration phase that uses the available alcohol challenge calibration data {ũ,ỹ} to optimally choose r 1 and r 2 . This takes the form of a second calibration phase optimization wherein we choose nonnegative and which minimize: where ũ*(t,r 1 , r 2 ) denotes the minimizer of equation (2.15) with ỹ = ỹ and regularization weights r 1 and r 2 . Minimizing J C,R chooses and in such a way that the fidelity of both the forward convolution and the deconvolution are taken into account.

A Nonlinear Pharmacokinetic Model for Alcohol Absorption, Metabolism, and Elimination
Ethanol, the form of alcohol contained in alcoholic beverages, ingested into the body is absorbed into the blood stream through the stomach and intestines. Ethanol, being highly miscible in water, and the human body, being composed of 70% water, ingested ethanol then ends up distributed throughout the body. Aside from the relatively small amount of ethanol processed out of the body by the kidneys and excreted through perspiration, the majority is metabolized into other waste products in the liver with the aid of a group of enzymes known collectively as alcohol and aldehyde dehydrogenase (ADH and ALDH) [8,21,44]. As is the case with most enzyme catalyzed reactions, at low substrate concentrations, the reaction exhibits first-order kinetics. At higher concentrations, however, the ability of ADH to catalyze the reaction becomes overwhelmed or saturated (this is essentially the basis for intoxication) and the kinetics become zero-order.
There are a number of models already in use that capture only the zero-order kinetics. That is because they are primarily used to determine if a subject is within the legal intoxication limit for operating a motor vehicle. As a result they are only accurate for relatively high BAC levels. The parameters in these models are determined empirically based on large scale studies using typically stratified (by gender, body mass, drinking history, etc.) population data. In our work we are interested in all levels of BAC including relatively low levels, especially for the purpose of calibrating our TAS models. In order to capture both the zeroand first-order metabolism/elimination kinetics, we model the absorption of ethanol through the gut using a linear term and the metabolism/elimination using a Michaelis-Menten [13] term that will capture the kinetics switching phenomenon. If we assume that initially there is no ethanol in the blood. Our model takes the form This optimization problem can be solved using a constrained (to maintain the non-negativity of the parameters) gradient projection method. The gradient of J can be computed once again using the adjoint method [24]. Indeed first note that, where .  for k = N i ,N i -1,... ,1,0, with z Ni = 2(u Ni -ũ Ni ). It follows that the gradient of J i , , can then be computed exactly with no truncation error, by first using (3.2) to compute , then using (3.3) to compute , and then evaluating We summarize the parameters that have to be fit in the following table.

Fitting the Nonlinear Pharmacokinetic Model
In practice, the parameters in the model given in (3.1) or (3.2) would be fit using population data (in the form of drinks consumed and time tagged measurements of BAC, or more likely, BrAC) from a large scale controlled study stratified by covariates such as body mass, sex, drinking history, etc. However, to demonstrate proof-of-concept, we present the results of fitting the model to a much smaller data set. One of the authors (SEL) served as the subject for this study (this was not considered to be human subjects research by the USC IRB). She 1) wore a WrisTAS™ 7 sensor over an 18-day period, during which she 2) collected breath measurements and 3) maintained a real-time drinking diary for all drinking episodes.
The WrisTAS™ 7 is worn like a wristwatch. Its sensing system directly measures the local ethanol vapor concentration over the skin surface at 5-minute intervals. The subject participated in an alcohol challenge consuming her first drink in the laboratory. BrAC was recorded every 15 minutes from the start of the drinking session until BrAC returned to 0.000. The subject then wore the sensor and consumed alcohol ad libitum for the following 17 days in the field. For each drinking episode, the subject would take BrAC readings every 30 minutes until the BrAC returned to .000. Figure 3 shows the entire 18 day TAC record along with the recorded BrAC measurements. The TAC readings are in units of mg/dl with the axis appearing on the left side of the plot. So as to be consistent with standard practice in the field, the BrAC measurements are reported in the units of % alcohol, or g/dl. The axis on the right hand side of the plot in Figure 3 is labeled in these units. The subject recorded 11 distinct drinking episodes which are noted in the figure.
The model (3.1) was sampled at 100 Hz. and fit to each of the 11 drinking episodes separately and then all together. The resulting values of the parameters are shown in Table 1 along with their means and standard deviations. We use the abbreviation Opt to denote the case of fitting all 11 episodes simultaneously. Note: To assess the quality of the fits, we then used the model given in (3.2) with the parameters given in Table 1 to simulate each of the episodes. For each of the 11 episodes, Table 2 shows the relative L 2 and L ∞ differences between the simulation and the data. It also shows the absolute differences in peak BrAC (P), time of peak BrAC (T) and area under the BrAC curve (A). To assess how well the model performed when it was used to simulate drinking episodes other than the one it was trained on, we simulated all 11 episodes using the model (3.2) with both the mean and Opt parameters. The results for the mean parameters are in the first five columns of Table 3 and  the results for the Opt parameters are in the last five columns of Table 3.
The mean differences indicate that the mean parameters seemed to perform somewhat better than the Opt parameters. Not surprisingly, neither did as well as when each episode's own training parameters were used. In general, the simulation was able to capture peak BrAC to within approximately 10% when the mean and Opt parameters were used, 5% when the training parameters were used. On average, in all three cases, it captured the time of peak to within about 20 minutes and the area under the BAC or BrAC curve to within approximately 5-10%.

Calibrating the TAC Sensor Data Analysis Software
We used estimated BrAC data obtained using (3.2) together with TAC data for the first drinking episode in Figure 3 to calibrate the data analysis software for the TAS sensor. In (3.2) we used parameter values for K, M, and α from Table 1. Specifically we used the parameter values obtained when the model (3.2) was trained using the first episode, subsequently this will be referred to as Method 7, the means of the parameters when (3.2) was trained on each of the 11 episodes individually (Method 8), and the parameters denoted in Table 1 by Opt (Method 9), the parameter values obtained when (3.2) was trained using the data from all 11 drinking episodes simultaneously. We compared these results to the calibrations that resulted when the recorded BrAC from episode 1 was used (i.e. the alcohol challenge data) and when simulated BAC was obtained by one of 6 standard zero order kinetics only models (Methods 1 -6). These models are described and compared in [8,14,22,25,31,33,49,50]  In Table 4 we display the values of the parameters in the model (2.2) -(2.5) and the regularized deconvolution performance index (2.15) obtained by using the 10 different methods for generating calibration BrAC data. As one might expect method 7 produced parameters closest to those obtained with the alcohol challenge data since it used values of K, M, and α fit to the very same BrAC alcohol challenge data. In Table 5 we display the results of using the software calibrated using BrAC data obtained by the different methods to invert or deconvolve BrAC estimates for the ten remaining drinking episodes. We display the mean absolute differences in peak BrAC (P), time of peak BrAC (T) and area under the BrAC curve (A) between estimates obtained with BrAC challenge data calibration and calibration using BrAC data generated by one of the nine methods, and their standard deviations.
In Figure 4, for drinking episode 1 in Figure 3, we have plotted the breath analyzer measured BrAC along with the estimated BAC computed using the nine different methods described above. The BrAC data quite clearly shows the transition from zero-order to first-order kinetics at approximately 2 hours after the beginning of the episode. Not surprisingly, the Michaelis-Menten scheme with parameters obtained by fitting the model (3.2) to the episode 1 BrAC (method 7) approximates the BrAC data quite well. As would be expected methods 1-6 display only the zero-order kinetics and are not that accurate a low concentrations. They do however capture the zero-order kinetics well, especially considering that they use population parameters and have not been fit to the data. The Michaelis-Menten models (Methods 8 and 9) that use the average values of the parameters from the 11 episodes and the optimal parameters for all 11 episodes taken all at once display the characteristic transition from zero-order to first-order kinetics, but overall appear to be less accurate than the linear models. Method 9 appears to out-perform method 8.
In Figure 5 we have plotted the forward convolution kernels, k(·) = k(q; ·), from (2.12) corresponding to the values of q obtained by calibrating using BrAC data from the alcohol challenge and via simulation using the methods 7, 8, and 9. In Figure 6 we have plotted simulations of the first drinking episode using the model (3.2) with parameters from Table 1 obtained by training on episode 1, the means of the parameters obtained by training on each of the 11 episodes individually and the parameters obtained by training on all 11 episodes simultaneously (Opt).
In Table 5 we see that the estimated BrAC obtained by the software was in general more accurate when it was calibrated using BrAC data generated by methods 1 -6 than it was with methods 7 -9. In Figure 5 we see that the convolution kernel resulting from the alcohol challenge data based calibration is quite similar to the one resulting from the calibration using BrAC generated by (3.2) with the parameters resulting from training based on episode 1. The kernels based on BrAC data from (3.2) with the mean and Opt parameters are quite different, but similar to each other, in particular with regard to peak BrAC. Finally, in Figure  6 we note that although as expected the simulation based on the Episode 1 parameters is closest to the BrAC observations, the simulations based on the mean and Opt parameters, are still reasonably close.

Discussion and Concluding Remarks
In general, from Table 2 we see that the model fits given in Table 1 fit the data well and that the Michaelis-Menten based models do a good job of capturing both the zero order kinetics at high BAC/BrAC and the first order kinetics at lower concentrations. For each of the 11 drinking episodes identified in Figure 3, when comparing the fit model to the data, the mean relative L 2 error was 24% ± 7% while the mean relative L ∞ error was 21% ± 5%. On average the fit models were able to capture the peak BAC/BrAC to within 0.003 percent alcohol and the time of the peak to within approximately 17 minutes. On average they captured the area under the BAC/BrAC curve to within 0.01 percent alcohol -hours. The fact that in the case of each of the 11 drinking episodes, the Michaelis-Menten model fits the data very well, but as can be seen from Table 1, there is significant variance in the values of the parameters, indicates that there are other un-modeled dynamics at play in the body's ingestion, metabolism and elimination of ethanol. These might include stomach contents, time of day, type of drink, etc. This would be an area for future study. Also, because of the nature of the data we had available to us in this proof-of-concept study, being from only a single subject, we expect that there would also be variance in parameter values from subject to subject, even when other factors such as sex, body mass, what was drank, stomach contents, etc. are controlled. Putting the Michaelis-Menten models on par with the linear zero-order kinetics models would require a large scale study with many subjects stratified according to the kinds of covariates listed above. This would be another important area for future research. The results in Table 3 give some indication of the robustness of the fits and the models might do when population parameter values are used. When the average (over the 11 fits) values of the parameters are used, the mean (again over the 11 drinking episodes) relative L 2 error in BAC/BrAC was 47% ± 13% while the mean relative L ∞ error was 37% ± 20%. On average, the peak BAC/BrAC was captured to within 0.007 and the time of peak to within approximately 20 minutes. When the optimal (over all 11 drinking episodes simultaneously) values of the parameters are used, the mean (again over the 11 drinking episodes) relative L 2 error in BAC/BrAC was 59% ± 27% while the mean relative L ∞ error was 47% ± 23%. The error in peak BAC/BrAC was approximately 0.008 and the error in the time of peak was approximately 17 minutes.
Although we are primarily interested in the ability of the model given in (3.1) or (3.2) to accurately predict BAC based on drink diary input, it is of some interest and it is instructive to compare the values we obtained for the Michaelis-Menten parameters given in Table 1 when we fit the model to data, to values for these parameters appearing elsewhere in the literature. For example, in [48] the authors estimated the parameters Mand Kusing data obtained by dosing a single male volunteer (age 50, wt. 63.5 kg, ht. 179 cm) with 95% ethanol and then drawing and assaying his blood several times an hour over a 6 hour period. The study was repeated 5 times with different initial doses of alcohol, but the conditions under which they were carried out, were to the best of the researchers' and subject's abilities, kept constant (pre-study fasting, standardized meals, times between meal and dosing, etc.). The authors then fit Michaelis-Menten clearance or elimination models to each of the trials and then computed average values for the parameters. They found K = 39.87 mg/(dl h) and M = 16.56 mg/dl., with a high degree of variability from trial to trial. Now while there are a number of factors that would be expected to affect these values (e.g. gender, heavy drinker vs teetotaler, weight, height, age, type of dose, stomach contents, etc.), essentially all of which are significantly different in this study than in our own, these values are still an order of magnitude larger than the ones we obtained in our fits as shown in Table 1. On the other hand, in one study described in [16], the authors used values for K ranging from 22.30 to 70.64 mg/(dl h) and for M from 0.80 to 4.6 mg/dl. In another study those same authors used K = 47mg/(dl h) and M = 38 mg/dl for men and K = 48 mg/(dl h) and M = 40.5 mg/dl for women, and in a third study, they used values for K ranging from 54.26 to 100.00 mg/(dl h) and M = 4.5 mg/dl. However, they also indicated that you could obtain similar results by doubling K and multiplying M by approximately 6.5. The model used in [51] and [52] was a two compartment model, one compartment being essentially all of the blood and body water and the second being the liver. They use a Michaelis Menten clearance term in the equation for the liver. In order to compare the values for Mand K that they used in those studies to the ones we obtained in Table 1, their two compartments must be combined into a single compartment. It is not possible to do this after the fact in any precise way because of the nonlinearity of the Michaelis Menten term. However, since the only elimination term in their entire model is the Michaelis-Menten term in the equation for the liver compartment, it is possible to obtain a rough estimate for the values they used for Mand Kthat would be comparable to ours by simply replacing the volume of the liver with the subject of our study's volume of body water. In this way we obtain K = 21.0 mg/(dl h) and M = 0.46 mg/dl.
An alternative way to compare the values we obtained for the Michaelis Menten clearance parameters, K and M, with the values appearing elsewhere in the literature that will let us look at just a single value rather than two, is by considering the slope in a Lineweaver-Burk, or double reciprocal plot. Before the advent of high speed computing and sophisticated numerical methods for solving nonlinear least squares optimization problems, Michaelis-Menten models were typically fit to data by applying standard linear regression techniques to a Lineweaver-Burk plot of the data. If we let R denote the alcohol clearance rate and C the concentration of alcohol, then the Michaelis-Menten model takes the form or equivalently where λ = M/K and ρ = 1/K. For our fits, the average over the 11 drinking episodes was λ = 1.24 h, whereas for the studies mentioned above, values of λ from .01 to 2.44 were observed.
From Figure 4 we see again that the when fit directly to the data (Method 7), the Michaelis-Menten model does a very good job at capturing the zero and first order kinetics and the transition from one to the other. When population parameters are used, the optimal parameters (Method 9) appear to be preferable to the averaged parameters. Among the strictly zero order kinetics models, Methods 1,3,5 and 7 appear to capture the high concentration behavior well, whereas Method 6 appears to be more accurate at lower levels of BAC/BrAC. Tables 4 and 5 and Figures 5 and 6 reveal what happens when the estimated BAC/BrAC is used to calibrate our transdermal models. Table 4 shows that not surprisingly (based on Figure 4) Method 7 yields almost the same values for the transdermal parameters as does the alcohol challenge BrAC data. Methods 6 and 9 also yield parameter values that are close. Figure 5 shows how the variance in parameter values in Table 4, and the error in fitting the BrAC data for Episode 1 shown in Figure 4 between the Michaelis-Menten model fit directly to the data (Method 7) and the ones using population parameters (Methods 8 and 9) translate to variation in the convolution kernel and the convolution model itself. Then when the resulting calibrated (using just Episode 1 data) transdermal model is used to deconvolve BrAC from TAC across all 11 field episodes, comparing error in peak BrAC, time of peak and area under the curve, Table 4 appears to indicate that Methods 3 and 6 yielded the most accurate estimates. In general, the Michaleis-Menten models seemed to perform less well than the strictly zero order kinetics models, however, they did do better than when the alcohol challenge data was used to calibrate the transdermal model. Indeed, with the alcohol challenge data, we found that μ |ΔP| = 0.0114 with σ |ΔP| = 0.0109, and μ |ΔT| = 0.49 with σ |ΔT| = 0.30. These findings may be the result of the fact that the strictly zero order kinetics models, being fit to a large population data set are more robust and less susceptible to outliers across a variety of drinking episodes than are our three Michaelis-Menten models which were only trained on either a single episode or on this dataset alone. Finally, as might be expected, Figure 6 shows that when one looks at the deconvolution of BAC/BrAC from the TAC for drinking episode 1, Method 7 and the alcohol challenge data (i.e. the BrAC) yield very similar results (Recall Figure 4), the next best is Method 9 followed by Method 8.
In general we conclude that using estimated BAC/BrAC based on drinking data computed from physiological based models for the ingestion, metabolism and elimination of alcohol to calibrate our transdermal models is comparable to using BrAC that has been measured directly using a breath analyzer. We note that it should be possible to use these models (with population parameters) as part of a smart phone app that would make it convenient for the subject or patient to enter his drinking data for a single episode in the field. The app could email a report to the clinician or researcher which could then be automatically synchronized with the report from the TAC sensor and then used to calibrate the transport models that form the basis for the transdermal sensor data analysis system, thus eliminating the need for a laboratory or clinic alcohol challenge.  Schematic diagram of a calibration scheme for data analysis software based on drinking data rather than alcohol challenge data. A nonlinear pharmacokinetic Michaelis -Menten alcohol clearance model is used to estimate BAC from a patient supplied drinking diary. The resulting estimated BAC is then used with TAC data from the biosensor to calibrate the models used to deconvolve BAC from TAC in subsequent drinking episodes for which no drink diary is available.  Contemporaneous TAC and BrAC data for 11 drinking episodes recorded by a single subject engaging in naturalistic drinking in the field. A comparison of estimated BAC for drinking episode 1obtained using Methods 1 -9 based on drink diary input. The triangles are the breath analyzer measured BrAC. It is not surprising that Method 7 which uses the Michaelis -Menten alcohol clearance model with parameters fit based on drinking episode 1data is the most accurate. The Michaelis -Menten model with the mean (over all 11drinking episodes) parameters (method 8) also does reasonably well. As one would expect, the Michaelis -Menten models out -perform the 0order kinetics based models at low BAC.  Forward convolution kernels calibrated using TAC data from the first drinking episode and either BrAC data obtained via an alcohol challenge or by using estimated BAC based on drinking diary data computed using the Michaelis Menten alcohol clearance model with parameters as in methods 7, 8, or 9. TAC and BrAC data for the first drinking episode in Figure 3, along with estimated BAC or BrAC deconvolved from the TAC data from the first drinking episode using a forward model that was calibrated using either BrAC data obtained via an alcohol challenge or estimated BAC based on drinking diary data generated by our nonlinear pharmacokinetic model with a Michaelis Menten alcohol clearance term with parameters determined by (1) training using data from drinking episode 1 (solid line) (i.e. method 7), (2) averaging the parameters obtained by training on each of the 11 drinking episodes in Figure 3 Table 1 Parameters obtained by fitting (3.2) to each of the 11 drinking episodes in Figure 3. The parameter K is in units of mg/(dl h), M is in units of mg/dl and α is in units of mg/(dl std drinks h). The parameters in the row marked Opt are ones fit using the data from all 11 drinking episodes simultaneously.  Table 2 Comparison of simulation and data for each drinking episode shown in Figure 3. In addition to relative L 2 and L ∞ differences between the breath measured BAC (BrAC) and model estimated BAC, we compare peak BAC  Table 3 Simulation results using mean and Opt parameters from  Table 4 Calibration results in the form of parameter values (i.e. q 1 ,q 2 ) for the forward transdermal transport model and regularization weights (i.e. r 1 ,r 2 ) for the inversion or deconvolution process using measured BrAC (from an alcohol challenge) and estimated BAC obtained using 9 different methods.  Table 5 Inversion results using the 9 different BAC estimation methods. We compare average (over the 11 drinking episodes) absolute difference in peak BAC (units: % alcohol or g/dl), time of peak BAC (units: hours), and area under the BAC curve (units: g h/dl) when computed using forward and deconvolution models calibrated using BAC estimated using one of the 9 methods and forward and deconvolution models calibrated using BrAC obtained during an alcohol challenge session.