University of Southern California Dissertations and Theses

Estimation of random input to semi-linear abstract parabolic systems with application to quantitative description of drinking behavior based on trans-dermal alcohol concentration

(USC Thesis Other)

Estimation of random input to semi-linear abstract parabolic systems with application to quantitative description of drinking behavior based on trans-dermal alcohol concentration

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Estimation of Random Input to Semi-linear Abstract Parabolic Systems With Application to Quantitative Description of Drinking Behavior Based on Trans-dermal Alcohol Concentration Shuang Li Advisor: I.G. Rosen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Applied Mathematics) August 2019 Copyright 2019 Shuang Li Acknowledgement Foremost, I would like to express my special thanks of gratitude to my advisor Prof. I. Gary Rosen for his support of my Ph.D. study, research, and especially the writing of this thesis. In the past three years, I have been working together with Professor Rosen, and I have learned so much from him. I really appreciate for his time and eort to guide me through graduate school. Besides my advisor, I would like to thank the rest of my thesis committee: Prof. Susan E. Luczak and Prof. Chunming Wang for their encouragement and insightful comments. Last but not the least, I would like to thank my parents: Shuping Li and Junxiang Li, for giving me the courage to pursue a Ph.D. degree, for always reminding me the power of education, and for always telling me to be a better person. Thanks for all your encouragement! i Contents 1 Introduction 2 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 A Brief Summary of Some Previous Work . . . . . . . . . . . . . 8 1.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.6 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Basic Denitions and Preliminaries 13 3 A Model for Alcohol Transport in the Human Body 18 4 Results for Linear Parabolic Systems 23 4.1 Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 The Abstract Cauchy Problem with Unbounded Input . . . . . . 27 4.3 Abstract Linear Initial Boundary Value Problems with Input and Output on the Boundary of the Spatial Domain . . . . . . . . . . 30 5 Abstract Semilinear Initial-Boundary Value Problems with In- put on the Boundary of the Domain 36 6 Sampled Time Formulation 44 7 Finite Dimensional Approximation and Convergence 47 8 Parameter Estimation 50 9 Application to our Alcohol Model 53 ii 9.1 Estimating Blood or Breath Alcohol Concentration from Biosen- sor Measured Transdermal Alcohol Concentration . . . . . . . . . 57 9.2 The Adjoint Method . . . . . . . . . . . . . . . . . . . . . . . . . 60 10 Numerical Studies for the Deterministic Model 63 10.1 Deterministic Model . . . . . . . . . . . . . . . . . . . . . . . . . 63 10.2 Finite Dimensional Approximation . . . . . . . . . . . . . . . . . 65 10.3 Experimental Data Description . . . . . . . . . . . . . . . . . . . 67 10.4 Estimating the Model Parameters Based on One Drinking Episode 68 10.5 Model Parameters Using Multiple Drinking Episodes . . . . . . . 71 10.6 Estimating Model Parameters with Regularization . . . . . . . . 75 11 Estimation of Random Dynamical Systems 80 11.1 The Prohorov Metric on a Set of Probability Measures . . . . . . 80 11.2 Estimation of Random Discrete Time Dynamical Systems . . . . 83 11.3 Using the Prohorov Metric Framework to Non-Parametrically Es- timate the Underlying Probability Measure . . . . . . . . . . . . 90 11.4 Parametrically Estimating the Probability Density Function . . . 97 12 Estimation of Random Parameters in a Distributed Parameter Model for the Metabolism and Transport of Ethanol in the Human Body Based on Transdermal Alcohol Concentration 101 13 Quantitatively Characterizing and Estimating Drinking Pat- terns and Behavior 109 13.1 Interpretation of Drinking Parameters . . . . . . . . . . . . . . . 109 13.2 Data Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 13.2.1 The Stochastic Model . . . . . . . . . . . . . . . . . . . . 111 13.2.2 The Simulator and Sampled Data . . . . . . . . . . . . . 114 iii 13.3 The Parameter Estimation Scheme for Parametrically Estimating the Probability Density Function . . . . . . . . . . . . . . . . . . 116 13.4 The Parameter Estimation Scheme for Estimating the Probabil- ity Measure Based on the Prohorov Metric Framework . . . . . 120 13.5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 121 13.5.1 Numerical Results for Parametrically Estimating the Prob- ability Density Function . . . . . . . . . . . . . . . . . . . 121 13.5.2 Numerical Experiment 1 . . . . . . . . . . . . . . . . . . . 122 13.5.3 Numerical Experiment 2 . . . . . . . . . . . . . . . . . . . 122 13.5.4 Numerical Experiment 3 . . . . . . . . . . . . . . . . . . . 123 13.5.5 Numerical Experiment 4 . . . . . . . . . . . . . . . . . . . 123 13.5.6 Numerical Experiment 5 . . . . . . . . . . . . . . . . . . . 124 13.5.7 Numerical Experiment 6 . . . . . . . . . . . . . . . . . . . 124 13.6 Numerical Results for Non-Parametrically Estimating the Drink- ing Parameters using Prohorov Metric Framework . . . . . . . . 125 13.6.1 Numerical Experiment 7 . . . . . . . . . . . . . . . . . . . 125 13.6.2 Numerical Experiment 8 . . . . . . . . . . . . . . . . . . . 127 14 Summary, Concluding Remarks and Discussion of Future Re- search 129 14.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 14.2 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 130 14.3 Avenues of Possible Future Research . . . . . . . . . . . . . . . . 130 Bibliography 132 iv List of Tables 10.4.1 Values of the model parameters and various error values when calibrating on the 11 dierent drinking episodes. . . . . . . . . . 69 10.4.2 Summary of data presented in Table 10.4.1 . . . . . . . . . . . . 70 10.5.1 Various error values when calibrating on the 11 dierent drinking episodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 10.5.2 Summary of data presented in Table 10.5.1 . . . . . . . . . . . . 73 10.6.1 Values of the model parameters when calibrating on the 11 dif- ferent drinking episodes . . . . . . . . . . . . . . . . . . . . . . . 76 v List of Figures 1.1.1 Alcohol Biosensor Devices: Left panel: The Giner WrisTAS. Right panel: The AMS SCRAM. . . . . . . . . . . . . . . . . . . 3 1.3.1 Michaelis{Menten saturation curve . . . . . . . . . . . . . . . . 5 1.3.2 Observed data (Wilkinson et al., 1977) versus model-predicted blood ethanol curves after ingestion of four dierent doses of ethanol in adult white male subjects . . . . . . . . . . . . . . . . 6 10.3.1 Experimental Data: BrAC and WrisTAS TM 7 alcohol biosensor TAC measurements. . . . . . . . . . . . . . . . . . . . . . . . . . 67 10.4.1 Estimated and Experimental TAC and BrAC using parameters estimated from Episode 1 . . . . . . . . . . . . . . . . . . . . . . 70 10.4.2 Estimated and Experimental TAC and BrAC using parameters estimated from Episode 4 . . . . . . . . . . . . . . . . . . . . . . 71 10.5.1 Estimated and experimental TAC and BrAC for drinking episode 1 using parameters estimated using drinking episodes 1-7 . . . . 73 10.5.2 Estimated and experimental TAC and BrAC for drinking episode 4 using parameters estimated using drinking episodes 1-7 . . . . 74 10.5.3 Estimated and experimental TAC and BrAC for drinking episode 5 using parameters estimated using drinking episodes 1-7 . . . . 74 10.5.4 Estimated and experimental TAC and BrAC for drinking episode 11 using parameters estimated using drinking episodes 1-7 . . . 75 10.6.1 Estimated and experimental TAC and BrAC using parameters estimated based on drinking episode 1 data with regularization. 77 10.6.2 Estimated and experimental TAC and BrAC using parameters estimated based on drinking episode 7 data with regularization. 77 vi 10.6.3 Estimated and experimental TAC and BrAC for drinking episode 1 using parameters estimated using drinking episodes 1-7 but now with regularization . . . . . . . . . . . . . . . . . . . . . . . 78 10.6.4 Estimated and experimental TAC and BrAC for drinking episode 5 using parameters estimated using drinking episodes 1-7 but now with regularization . . . . . . . . . . . . . . . . . . . . . . . 79 10.6.5 Estimated and experimental TAC and BrAC for drinking episode 11 using parameters estimated using drinking episodes 1-7 but now with regularization . . . . . . . . . . . . . . . . . . . . . . . 79 13.2.1 Simulated BrAC with density parameters = 5; = 1; = 1 . . 116 13.2.2 Simulated TAC with density parameters = 5; = 1; = 1 . . 116 13.6.1 Left:Approximation of Distribution ofd using Prohorov Metric; Right: Actual Distribution ofd used to Generate Data . . . . . 126 13.6.2 Left:Approximation of Distribution ofd using Prohorov Metric; Right: Actual Distribution ofd used to Generate Data . . . . . 128 vii Abstract We developed a deterministic semilinear parabolic system to describe alco- hol metabolism and transport in the human body based on transdermal alcohol concentration. The novelty of this model is that we used an ordinary dierential equation to describe the elimination process of blood alcohol and a partial dier- ential equation to describe the transport of alcohol molecules through the skin. In this model, we assume that ingested alcohol is immediately absorbed into the blood and the amount of alcohol excreted through the skin is proportional to the amount of alcohol in the blood. We also built a computational scheme to train individual's parameters in the alcohol metabolism and transport model from one or multiple drinking episodes of transdermal alcohol concentration data. With trained parameters, we can evaluate blood alcohol concentration from transdermal alcohol concentration. Later in this thesis, we constructed a random discrete time dynamical system by assuming input parameters in the deterministic semilinear parabolic system are random. Therefore, the drinking behavior estimation problem is formulated as a problem of identifying random parameters in such system. By establishing a Prohorov metric framework to non-parametrically estimate the distribution of the random parameters and a parametric framework to estimate the distribution of the random parameters, we provided a possible way to quantitatively study people's drinking behaviour from transdermal alcohol concentration collected from non-intrusive wearable biosensor. 1 Chapter 1 Introduction 1.1 Background It is widely known that over-drinking can cause various health problems. A large amount of research concerned with drinking behaviors has been done in recent decades. Although drinking volume plays an important role in drinking behaviors, drinking pattern has been linked to drinking behaviors as well. Research on drinking patterns is important, because the way alcohol is con- sumed has dierent eects given the same volume (e.g. the dierent impact of drinking 1 drink every day for 7 days vs. 7 drinks in 1 day). Drinking pattern is generally classied as one of the following: occasional binge drinking, heavy social drinking, alcohol abuse and so on. However, the data for this type of study is usually hard to collect and un- reliable. First of all, it's almost impossible to track people's drinking pattern. In most similar research, respondents are asked to report their usual drinking frequency and quantity of consumption. However, the reliability of this type of data is problematical, since there is no simple method to measure how much ethanol is in dierent drinks. Secondly, records of drinking volume and fre- quency after respondents are intoxicated are unreliable. A more precise and reliable method to collect data is blood alcohol con- centration (BAC) which indicates drinking consumption instantly. The most commonly used device to estimate BAC is the Breathalyzer, but it requires a small breath sample for every estimation. So, it is not practical to measure BAC at short intervals and the data collection process will likely have an eect 2 on drinking pattern. Therefore, using a trans-dermal alcohol sensor (for example, in the form of a wearable bracelet) would be a non-intrusive and passive way to obtain alcohol consumption data, that could be used to monitor BAC over relatively long periods of time. Several devices, including the WrisTAS TM (Giner, Inc., Newton, MA) and SCRAM c (Alcohol Monitoring Systems, Inc., Denver, CO) have been developed to measure transdermal alcohol concentration (TAC) at essentially continuous rates for weeks at a time (See Figure 1.1.1 below). Figure 1.1.1: Alcohol Biosensor Devices: Left panel: The Giner WrisTAS. Right panel: The AMS SCRAM. 1.2 Purpose Our purpose is to make use of TAC data collected from a trans-dermal alco- hol sensor to study an individual's drinking pattern quantitatively, to establish a mathematical model to describe the movement of ethanol molecules in the human body and to solve for the parameters in the mathematical model. However, deterministic models are not sucient to study drinking pattern quantitatively. Firstly, there is a need for many parameters to be determined in a deterministic well-structured model. Usually, determination of parameters requires a large number of experiments. On the other hand, in order to establish a pattern, a stochastic model can describe the distribution by using very few 3 parameters. Secondly, models concerning alcohol metabolism with ingestion of a single drink are not feasible, because ingestion of multiple drinks in a short span of time is common. Thirdly, a lot of factors that aect people's drinking patterns are dicult or impossible to model, such as emotion, occasion and event. Therefore, stochastic models are a better t to study drinking pattern compared to deterministic models. 1.3 Method To accomplish this study, we will use a partial dierential equation(PDE) and ordinary dierential equation(ODE) hybrid system to describe the move- ment of ethanol molecules in the human body. Wang(2006) has established a model of alcohol digestion with multiple compartments which is also a hybrid system. The advantage of our model is that it connects BAC and TAC via a boundary condition of the PDE, instead of adding BAC to the PDE as in the previous model. Furthermore, in order to study drinking patterns, we con- sider a stochastic input to the PDE, which forms a stochastic partial dierential equation(SPDE). In order to formulate the dynamics of the alcohol concentration in the blood as an ODE, we will brie y explain the biochemical reactions involved in alco- hol metabolism. Mainly, there are two enzymes involved in the process of al- cohol metabolism, alcohol dehydrogenase (ADH) and aldehyde dehydrogenase (ALDH). First, ADH oxidizes ethanol into acetaldehyde. Then acetaldehyde is metabolized by ALDH to acetate, which can be processed out of the body. Therefore, the rate of alcohol metabolism will be controlled by the concentra- tion of enzymes. And the most commonly used model to describe this type of 4 biochemical reaction is so-called Michaelis-Menten kinetics given by, v = dw dt = V max w K M +w (1.3.1) wherev is the rate of alcohol metabolism, w is the alcohol concentration in the blood, V max represents the maximum rate of alcohol metabolism achieved at the saturation concentration of the enzymes, and the Michaelis constant K M is the concentration of enzymes such that the reaction rate is half of V max . Figure 1.3.1: Michaelis{Menten saturation curve Figure 1 shows how Michaelis-Menten kinetics aect the blood ethanol con- centration after single dose alcohol ingestion, where the vertical axis is the absolute value of reaction ratejvj and the horizontal axis is alcohol concen- tration w. When the concentration w is small at the beginning, we nd that jvj = Vmaxw K M from equation (1.3.1), so the absolute value of reaction rate grows approximately linearly with respect to w. As the concentration increases, we obtainjvj = V max , which implies that the absolute value of the reaction rate is approaching the maximum rate V max . Equation(1.3.1) indicates an expo- nential decrease for ethanol concentration when ethanol concentration is much lower than the saturation concentration. However, as the ethanol concentration approaches the saturation concentration, ethanol concentration will decrease linearly, i.e. at an approximately constant rate. 5 Figure 1.3.2: Observed data (Wilkinson et al., 1977) versus model-predicted blood ethanol curves after ingestion of four dierent doses of ethanol in adult white male subjects Alcohol transport through the skin is modeled by a diusion equation, be- cause we assume that the movement of ethanol molecules is a diusion process. The nonlinear ODE describing the reaction and the linear PDE describing dif- fusion are linked through the boundary condition where the ux of alcohol at the skin boundary is proportional to the alcohol concentration in the blood. Therefore, we have a hybrid system, which can be generalized to an abstract system with in general unbounded input and output. After deriving the model and its well-posedness properties, we develop an abstract approximation framework and convergence theory for formulating and solving this type of estimation problem. Our approach relies heavily on two papers: (1) Banks's framework for the approximation of nonlinear functional dierential equation control systems, which shows that under certain assump- tions on the nonlinearity, the existence and uniqueness of the solution to the innite dimensional systems can be proved. And (2) Curtain and Salamon's 6 theory for constructing nite dimensional compensators for a class of innite dimensional systems with unbounded input operators. While our eort here is similar in spirit from the Curtain and Salmon treatment, it is dierent in that we have not only unbounded input but also nonlinearities in the model. Because of the uncertainties in how a subject ingests alcohol, as we noted above, a stochastic model for the input to the system is better suited than a deterministic model. Since ingestion of alcohol will almost immediately cause a change in ethanol concentration in the body water, it is reasonable to assume that alcohol concentration in the blood will be increased by an amount pro- portional to the amount of alcohol ingested. For example, one possible model for the input is to assume that alcohol is consumed at random magnitudes fA i ;i = 1; 2;:::g at the jump times of independent Poisson processes N(t). Such an input might take the form u(t) = 1 i=1 A i (!)1 t=i(!) (1.3.2) In order to estimate the parameters in a stochastic input such as this, rather than estimate the input directly, we formulate the problem of characterizing and identifying drinking patterns as the one of estimating parameters associated with random variables appearing in the model. In this way we did not have to develop a theory for estimating random inputs, and instead could rely on, and appropriately modify, existing treatments for estimating random parameters in abstract dynamical systems. One approach is based on the Prohorov metric on the space of probability measures (see, for example, [2, 5]) and the other involves the estimation of the probability distribution in the form of a probability density function for random parameters in abstract parabolic PDEs (see, for example, [29, 30, 31, 32]). We also consider two dierent statistical models, the so-called naive pooled data model (see, for example, [2, 5])) and a mixed eects 7 model (see, for example, [15, 16, 17]). Mixed eects models take into account uncertainty both within and across subjects. 1.4 A Brief Summary of Some Previous Work It has been known since the 1930's that the concentration of ethanol in sweat and perspiration is positively correlated with the concentration of ethanol in all of the body water and in particular with the blood (see, for example, [23]). More recently, in the past two or three decades, there has been expanded re- search activity on this topic (see, for example, [34, 35, 36, 37]). This research has been driven 1) by the development of new transdermal technology that al- lows one to quantitatively measure the amount of ethanol passing through and exiting the surface of the skin, for example, the WrisTAS 7 Transdermal Alco- hol Sensor by Giner, Inc. in Waltham Massachusetts, 2) the increasing need to monitor people with alcohol use disorder, for example using the Secure Con- tinuous Remote Alcohol Monitoring System (SCRAM) by Alcohol Monitoring Systems (AMS) in Littleton, Colorado, and 3) the demand for wearable health monitoring technology by the general public as evidenced by devices being de- veloped by companies such as Bactrac in San Francisco, California, and Milo Sensors in Santa Barbara, California. The underlying chemical principle behind the majority of these devices is an oxidation{reduction (or redox) reaction that yields 4 electrons for each molecule of ethanol. This produces a current which can then be measured. Alternatively, some of the newer devices employ an enzymatic reaction. Regardless of the un- derlying technology, however, the problem is, as it is with virtually all biosen- sor devices, to then reliably, accurately, consistently, and continuously convert whatever analog signal the sensor measures into the familiar and commonly used measure of intoxication in the human body, blood alcohol concentration 8 (BAC) or breath alcohol concentration (BrAC). Until these new transdermal devices came along, the level of alcohol in the body was measured using the breath analyzer. The technology in the breath analyzer is based on Henry's Law from elementary chemistry which is used to describe the exchange of gases between circulating pulmonary blood and alveolar air (see, for example, [21]). The breathalyzer is considered to be reasonably robust across the general pop- ulation, enough so that it can serve as the basis for arrest for driving under the in uence (DUI) of alcohol. Unfortunately the transport of ethanol through the skin and then its mea- surement by an analog sensor is both physiologically and technology more chal- lenging than simply measuring the amount of ethanol being exhaled from the lungs. Even more signicant is the fact that the alcohol signal being excreted transdermally or through the skin has been observed to be much more directly aected by several factors that vary between subjects (e.g., skin layer thickness, porosity and tortuosity, etc.) and even between individual drinking episodes by the same subject (e.g., body and ambient temperature, skin hydration, va- sodilation). Consequently the challenge in making transdermal alcohol sensing devices feasible for research, clinical, and even recreational applications is the development of theory and algorithms that reliably and consistently convert the biosensor measured transdermal alcohol concentration (TAC) to either BAC or BrAC. In earlier work by Rosen and Luczak's research group at the University of Southern California (USC) ([13, 18, 27]) deterministic nonlinear least squares techniques have been used to t rst principles physics and physiologically based distributed parameter or partial dierential equation models for the forward transport of ethanol from the blood through the skin and its eventual mea- surement by the transdermal sensing device. These studies have included the 9 formulation of the models, the development of nite dimensional approxima- tion schemes and convergence arguments, gradient based optimization methods using the adjoint method [22] to eciently and accurately compute gradients, and extensive numerical studies involving both simulated and actual clinical data to demonstrate the ecacy of their approach. The models they developed primarily used either parabolic or hyperbolic PDEs to describe the diusion of ethanol through the interstitial uid between cells in the epidermal layer of the skin. The models also involved unbounded input and output. Convergence arguments for the nite dimensional approximations were based on the Trotter Kato theorem from linear semigroup theory [20, 25]. The group at USC then developed inversion or blind deconvolution methods that use the t model to estimate BAC or BrAC (which one depends upon which type of data was used to t or train the model) from the biosensor measured TAC. The USC group next turned their attention to developing models that did not have to be calibrated to each individual subject. They t population models by considering the parameters in the PDEs to be random variables and then estimating their distributions. This work is described in the papers ([29, 30, 31, 32]) and is based on ideas in the two papers [19, 28]. In these treatments, the authors eectively consider the random parameters as if they were additional space variables in the PDE and then look at the weak form of the equation in appropriately constructed Bochner spaces with ranges in the original Gelfand triple of spaces in which the PDE is set. These Bochner spaces themselves form a Gelfand triple and the weak form is shown to be well-posed in these spaces. Generation and approximation results from linear semigroup theory are then used to establish a convergence and consistency result for the estimated probability density functions in the spirit of the ones in [2, 5]. The researchers were also able to extend their earlier deconvolution results to these new systems 10 and use the distribution of the random parameters to construct credible bands for the estimated BAC/BrAC signal and credible intervals for a number of statistics of interest to alcohol researchers including peak BAC/BrAC, time of peak BAC/BrAC, and area under the BAC/BrAC curve. 1.5 Outline An outline of the remainder of the thesis is as follows. In Chapter 2, we rst review some of the basic denitions and preliminary results that will be of use to us later in the thesis. We then in Chapter 3 establish a model in the form of an abstract semilinear hybrid ODE and PDE based parabolic system with input and output on the boundary that describes how alcohol is metab- olized and transported in the body. In Chapter 4 we review existing results for abstract linear parabolic systems including for those with unbounded input and in particular input on the boundary of the spatial domain. In Chapter 5, we extend the results in Chapter 4 to the semilinear case. In Chapter 6 we discuss the sampled time version of the model and in Chapter 7 we consider nite dimensional approximation and convergence. Chapter 8 contains our re- sults applying the theory developed in the earlier chapters to the deterministic parameter estimation problem. In Chapter 9 we apply the theory developed in Chapters 5 through 8 to the alcohol metabolism and transport model derived in Chapter 3. Chapter 10 has the results and a discussion of our numerical studies for the deterministic parameter estimation problem. In Chapter 11 we review some recent results for the estimation of random dynamical systems and in Chapter 12 we apply these ideas to the estimation of drinking patterns based on the model derived earlier in the thesis. Chapter 13 has our numerical studies for the stochastic estimation problem and we conclude the thesis in Chapter 14 with a summary and a discussion of some future avenues of research. 11 1.6 Notation Here, we present the largely standard notation used in this thesis. For a closed interval I in R and a Banach space X, L p (I;X) denotes the Banach space consisting of equivalence classes of strongly Lebesgue-measurable func- tions f : I ! X with R I jjfjj p X < 1. The induced norm on L p norm is jjfjj Lp = ( R I jjfjj p X ) 1=p when 1 p < 1. We denote by C(I;X) the Ba- nach space of continuous functions from I into X with the usual supremum norm. We letC k (I;X) denote the space of all continuous functions fromI toX whose kth derivatives are also continuous. The above notations are shortened to L p (I), C(I) and C k (I) respectively when X = R n . For Banach spaces X and Y ,L (X;Y ) denotes the space of all bounded linear operators from X to Y . In the case that X = Y , above notation can be shortened toL (X). For the most part, we use X and H to represent Banach space and Hilbert space respectively. For a linear operatorA, we employ(A);(A);(A) to denote the spectrum, point spectrum and resolvent set in the complex plane C of A, and D(A) andR(A) denote its domain and range. 12 Chapter 2 Basic Denitions and Preliminaries In this section, we will summarize some of the basic denitions and theorems we use in the following sections, so that we can refer to this section when we need to introduce these ideas and results in what is to follow. Unless otherwise specied, throughout this section, we let X denote a real or complex Banach space andfS(t) : 0t<1g be a family of operators inL (X) parameterized by a nonnegative real number t. Denition 2.0.1. (Gelfand Triple) Let H be a Hilbert space and let V be a dense subspace ofH such thatV can be endowed with a topological vector space structure for which the inclusion mapi is continuous, or equivalently, bounded. Identifying H with its dual space H via the Riesz map, the adjoint to i is the map i :H =H !V : It follows thatH is densely and continuously embedded in V and the duality pairing betweenV andV , denoted by<;> V ;V is then compatible with the inner product on H, in the sense that <u;v> V ;V =<u;v> H wheneveru2H =H V V andv2V H. The tripleV HV is then referred to as a Gelfand triple. Denition 2.0.2. (Sesquilinear Form)A sesquilinear form on a vector space H 13 over a eld F is a map a(;) :HH!F which is conjugate linear in the rst argument, and linear in the second, i.e. a(v 1 ;cw 1 ) =ca(v 1 ;w 1 ); a(v 1 ;w 1 +w 2 ) =a(v 1 ;w 1 ) +a(v 1 ;w 2 ) a(cv 1 ;w 1 >= ca(v 1 ;w 1 ); a(v 1 +v 2 ;w 1 ) =a(v 1 ;w 1 ) +a(v 2 ;w 1 ) for all c2F and v 1 ;v 2 ;w 1 ;w 2 2H. Denition 2.0.3. (Semigroup) A family of bounded linear operatorsfS(t) : t 0g on X is called a semigroup on X if (i)S(0) =I, (ii)S(t +s) =S(t)S(s) for every t;s 0. Denition 2.0.4. (Innitesimal Generator of Semigroup) LetfS(t) : t 0g be a semigoup on X, a linear operator A dened by Ax = lim t!0 + S(t)xx t with its domain of denition D(A) :=fx2X : lim t!0 + S(t)xx t existsg is called the innitesimal generator of the semigroupfS(t);t 0g. Denition 2.0.5. (C 0 -Semigroup) A semigroup fS(t) : t 0g is called a strongly continuous or C 0 -semigroup on X if lim t!0 + jjS(t)xxjj X = 0; for x2X Denition 2.0.6. (Analytic Semigroup) LetfS(t);t 0g be a C 0 -semigroup 14 on X with innitesimal generator A. fS(t);t 0g is said to be an analytic semigroup if (i)for some 0<<=2, the bounded linear operator S(t) :X!X can be extended to t2 , =f0g[ft2C :jarg(t)j<g and the usual semigroup conditions hold for s;t2 : S(0) = I;S(t +s) = S(t)S(s), and for each x2X, S(t)x is continuous in t; (ii)For all t2 nf0g, S(t) is analytic in t in the sense of the uniform operator topology. Denition 2.0.7. (Resolvent Set and Resolvent Operator) The resolvent set (A) of the operator A is dened as (A) =f2C :R (A) = (IA) 1 exists and boundedg and R (A) is called resolvent operator associated with A. The Hille-Yosida theorem provides necessary and sucient conditions for an operator to be the innitesimal generator of a strongly continuous semigroup. Theorem 2.0.1. (Hille-Yosida) LetA be a linear operator dened on a subspace D(A) of a Banach spaceX,M > 0 and!2R, thenA generates aC 0 semigroup fS(t);t 0g that satisesjS(t)jMe !t if and only if (i)A is closed and densely dened in X, i.e. D(A) =X, (ii)For real >!, we have 2(A) and R (A) satises jR (A) n j M (!) n ; n = 1; 2; Theorem 2.0.2. (Hilbert Projection Theorem) Let H be a Hilbert space and 15 M a closed subspace of H, there exists a unique point y2M for whichkxyk is minimized over M. The Trotter-Kato Theorem is a fundamental result on the approximation of semigroups. Theorem 2.0.3. (Trotter-Kato) Let H and H N be Hilbert spaces such that H N H. Let P N : H ! H N be an orthogonal projection of H onto H N . Assume P N x ! x as N ! 1 for all x 2 H. Let A N ;A be innitesimal generators of C 0 semigroups S N (t);S(t) on H N ;H respectively satisfying (i)There exists M;! such thatjS N (t)jMe !t for each N, (ii)There exists D dense in H such that for some , (IA)D is dense in H and A N P N x!Ax for all x2D. ThenS N (t)x!S(t)x uniformly int on compact intervals [0;T ] for eachx2H. Denition 2.0.8. A closed densely dened operator A is said to be of type (!;M) if there exist 0<! < and M 1 such thatf :jargj>!g(A) andjj(A) 1 jj M for < 0 and if there exists a number M such that jj(A) 1 jjM holds injargj>! + for all > 0. Theorem 2.0.4. (Necessary and Sucient Condition for the Generation of an Analytic Semigroup) IfA + is of type (!;M) for some real , where 0 ! < =2, then for each > ,f : Re > g (A) and in this half planejjA) 1 jj C(1 +j j) 1 , for some C > 0. Conversely if A is a densely dened closed operator satisfying the above bound then there exists a with < and !2 [0;=2] such thatA + is of type (!;M). IfA + is of type (!;M) where !2 [0;=2] and and M 1, then A is the generator of an analytic semigroup. Theorem 2.0.5. (Young's Inequality) If a and b are nonnegative real numbers 16 and p and q are real numbers greater than 1 such that 1 p + 1 q = 1, then ab a p p + b q q : Theorem 2.0.6. (Cauchy-Schwarz Inequality) For all vectors u and v of an inner product space it is true that j< u; v>j 2 < u; u>< v; v> where <;> is the inner product. Theorem 2.0.7. (Schmidt Inequality) If p n (x) is a polynomial of degree n = 1; 2; 3, then Z b a [Dp n (x)] 2 dx 4k n (ba) 2 Z b a [p n (x)] 2 dx; where k 1 = 3;k 2 = 15 and k 3 = 1 2 (45 + p 1605) 42:6. 17 Chapter 3 A Model for Alcohol Transport in the Human Body As discussed in [Kovatchev, et al., [7]], ethanol pharmacokinetics models have been developed from classic models including the Widmark Model [Wid- mark, [24]] to an enzyme-enhanced chemical reaction model with limited supply (Umulis et al., [39]). It is now widely accepted that the alcohol elimination in the liver follows the zero order and rst order Michaelis-Menten kinetics. The ethanol elimination process is governed by two main hepatic enzymes known as alcohol dehydrogenase (ADH) and cytochrome (CYP2E1)(AW, [1], J-W et al., [40]). Therefore, the rate of decline in the alcohol elimination in the liver with time is nonlinear and depends on the theoretical maximal decay rate V and the Michaelis-Menten constant K which indicates the alcohol concentration at V=2. So if we denote alcohol ingestion by u(t) and blood alcohol concentration by w(t), we have the following equation describing blood alcohol concentration with alcohol ingestion described by the input signal u, dw(t) dt = Vw(t) K +w(t) +bu(t);t> 0 Roughly 95% of the alcohol will be eliminated by oxidation in the liver, and the rest will be excreted through breath, sweat and urine. Specically, a portion of the ethanol molecules in the blood vessels will be transported through the boundary between skin and blood vessels. Furthermore, the spread of ethanol molecules in the skin can be considered as a diusion process. Thus, by denoting alcohol concentration at depth , 0 l in the skin(skin is assumed to be l cm thick) and time t, t > 0 as x(t;), we get the following hybrid system 18 with ODE and PDE to describe alcohol pharmacokinetics including ingestion, elimination and distribution: @x @t (t;) =D @ 2 x @ 2 (t;); 0l;t> 0 dw(t) dt = Vw(t) K +w(t) +bu(t);t> 0 D @x @ (t; 0) =x(t; 0) D @x @ (t;l) =w(t) x(0;) = 0; 0l w(0) = 0 y(t) =rx(t; 0) (3.0.1) In order to relax some of the restrictions on the solutions, we will study the weak formulation of the above system. Let H k ([0;l]) be Sobolev spaces W k;2 ([0;l]) which are sets of all functions f on [0;l] such that the functionf and its weak derivatives up to orderk are inL 2 ([0;l]) and have a niteL 2 norm when k = 1; 2. DeneH :=L 2 ([0;l])R,V :=H 1 ([0;l])R andW :=H 2 ([0;l])R, and their inner products as < ^ x; ^ > H = R l 0 x d +wr, < ^ x; ^ > V = R l 0 x d + R l 0 x 0 0 d +wr and < ^ x; ^ > W = R l 0 x d + R l 0 x 0 0 d + R l 0 x 00 00 d +wr, re- spectively, where ^ x = (x;w) and ^ = ( ;r). Denote H 1 0 ([0;l]) := f j 2 H 1 ([0;l]); (0) = 0; (l) = 0g and H 1 L ([0;l]) :=f j 2 H 1 ([0;l]); (0) = 0g. The inner product <; > without specifying the space will be the L 2 inner product<x; >= R l 0 x d, we can write the weak formulation of system (3.0.1) as: < _ ^ x; ^ > V ;V +a(^ x; ^ ;q) =f(w;u;q)r; for all ^ 2V (3.0.2) 19 where a(^ x; ^ ;q) =D<x 0 ; 0 >w (1) +x(0) (0) f(w;u;q) = kw M +w +bu and <; > V ;V now denotes the natural extension of the H innerproduct to the duality pairing betweenV andV and we have suppressed showing explicit dependence on t wherever possible. Denition 3.0.1. (Weak Solution) A weak solution to system(3.0.1) is a vector function ^ x(t)2C([0;T ];V )\W 1;2 ([0;T ];V ) such that system(3.0.2) holds for every t. Denition 3.0.2. (Classical Solution) A classical solution to system(3.0.1) is a vector function (x(t;);w(t)), where x(t;) dened on D = [0;T ] [0;l] such thatx(t;)2C 2 (D) andw(t)2C 1 ([0;T ]) and moreover system(3.0.1) holds at every (t;)2D. Theorem 3.0.1. The system(3.0.1) is equivalent to system(3.0.2) in the sense that every classical solution to system(3.0.1) is also a weak solution to sys- tem(3.0.2), and every weak solution to system(3.0.2) which satises the regu- larity conditions in Denition 3.0.2 is also a classical solution to system(3.0.1). Proof. Suppose ^ x = (x;w) is a solution to the system(3.0.1) and ( ;r)2V . By applying integration by parts, < ( _ x; _ w); ( ;r)> V ;V =D Z l 0 @ 2 x @ 2 d + ( kw M +w +bu)r =D Z l 0 @x @ 0 d +D @x @ l 0 kw M +w r +bur =D Z l 0 x 0 0 d +w (l)x(0) (0) kw M +w r +bur 20 which is equation (3.0.2). Conversely, suppose (x(t;);w(t)) is a solution to system(3.0.2) andx(t;)2 C 2 (D) and w(t)2C 1 ([0;T ]), < ( _ x; _ w); ( ;r)> V ;V =D< @ 2 x @ 2 ; >D @x @ l 0 +w (l)x(0) (0) + ( kw M +w +bu)r =D< @ 2 x @ 2 ; >D @x @ (l) (l) +D @x @ (0) (0) +w (l)x(0) (0) + ( kw M +w +bu)r Let 2H 1 0 ([0;l]) and r = 0, from equation above we can get < ( _ x; _ w); ( ;r)> V ;V =D< @ 2 x @ 2 ; > V ;V ; or < _ xD @ 2 x @ 2 ; > V ;V = 0 Since H 1 0 ([0;l]) is dense in L 2 ([0;l]), we obtain _ x =D @ 2 x @ 2 Now let 2H 1 L ([0;l]) and r = 0, we have D @x @ (0) H H H (0) =x(0) H H H (0) Let 2H 1 ([0;l]) and r = 0, we nd that D @x @ (l) H H (l) =w H H (l) 21 Plugging in all the above equations to equation (3.0.2), we nd that _ w A r = ( kw M +w +bu) A r Now, it is clear that system(3.0.1) is equivalent to system(3.0.2). Therefore, we can study a group of parabolic systems which generalizes system(3.0.2) in order to study system(3.0.1). In the next section, we will demonstrate how to nd the solution of this type of parabolic system. 22 Chapter 4 Results for Linear Parabolic Systems In this section, we will rst state the family of abstract semilinear parabolic systems that we will discuss in this thesis. We then review the results for the linear problem by specifying the requisite abstract operators and spaces. In sub- sequent chapters we will discuss 1) how to extend these results to the semilinear case, 2) the sampled time formulation, 3) nite dimensional approximation and convergence, 4) the application of these results to parameter estimation, and 5) the specic form these results take in the context of our alcohol metabolism and transport model. We will then discuss the results of our numerical studies before turning our attention to the problem of estimating stochastic input and drinking patterns. Of primary interest to us in this thesis is a class of abstract initial-boundary value problems. Unless specied explicitly, in this section, the following nota- tions are used. Let U = R and Y = R , W;V and H be Hilbert spaces such thatW;V ,!H with dense and continuous embeddings. TakingH as the pivot space, we can get that V ,!H ,!V and W ,!H ,!W are Gelfand triples. Let <; > denote the H standard inner product, andjj;jjjj denote norms on H and V respectively. Let Q be a compact metric space. Hence, system (3.0.1) in the last section can be generalized to the following 23 abstract initial-boundary value problem, _ '(t) = (q)'(t) +F (t;'(t);u(t);q);t> 0 (q)'(t) =G(t;'(t);u(t);q) '(0) =' 0 (q) y(t) =C(q)'(t) (4.0.1) where (q)2L (W;H), (q)2L (W;U) and C(q)2L (V;Y ), F (t;';u;q) : RWU!H and G(t;';u;q) :RWU!U. Before we consider the semilinear system(4.0.1) we review what has been done in the linear case. First, we discuss linear homogeneous Cauchy problems. Then we consider the nonhomogeneous problem with unbounded input in Sec- tion 4.2. Finally, we see how this theory can be used when the input is on the boundary of the spatial domain in Section 4.3. In the next chapter, we estab- lish analogous results for a semi-linear abstract initial-boundary value problems with unbounded input and output, i.e. system(4.0.1). Since we will discuss dif- ferent abstract problems in the following sections and chapters, we will specify the problem and the spaces where the problem is set at the beginning of each section . 4.1 Operators Since we need to deal with a number of dierent operators which act on dif- ferent spaces, we want to summarize all these operators and the spaces that they act on in this section. The main purpose is to reduce any confusion that simi- lar operators might cause and to demonstrate the relationship between similar operators. We start by dening A : V ! V based on a sesquilinear form a(;;q) : 24 VV !C that satises the following assumptions. Assumption. A1(Continuity) For q;e q2Q, we have for all '; 2V ja('; ;q)a('; ; ~ q)jd(q; ~ q)jj'jjjj jj where d(;) denotes any p-metric on R p . Assumption. A2(Coercivity) There exist c 1 > 0 and some real number 0 such that for q2Q, '2V , we have a(';';q) + 0 j'j 2 c 1 jj'jj 2 Assumption. A3(Boundedness) There exists c 2 > 0 such that for q 2 Q, '; 2V , we have ja('; ;q)jc 2 jj'jjjj jj Under these assumptions, a(;;q) denes a bounded linear operatorA(q) : V !V such that a('; ;q) =<A(q)'; > V ;V where '; 2V . Furthermore, it can be shown (Banks and Ito, [6], Tanabe, [38]) thatA(q) restricted on the set D(A(q)) =f'2V :A(q)'2Hg so that A(q) =A(q) onD(A(q)) is closed, is densely dened on H and is the innitesimal generator of an analytic semigroup S(t;q) on H. So we dene A(q) :D(A(q))!H 25 by restrictingD(A(q)) =f'2V :A(q)'2Hg. Therefore,A(q) is a restriction ofA(q). It can further be shown thatA(q) : V ! V is the innitesimal generator of an analytic semigroup onV that is an extension of the semigroupS(t;q) on H(see Tanabe, [38]). Banks and Ito [6] also shows that A(q) :D(A(q))H! H can further be restricted to an operator which is the innitesimal generator of an analytic semigroup on H. ~ A(q) is only used in Section 4.2 when we solve the general Cauchy problem and dened as follows. Since the operator A(q) is densely dened, closed and the innitesimal generator of a C 0 semigroup on H, it has an adjoint operator A(q) :D(A(q) ) H! H, which is densely dened and closed. Dene the spaceZ to be the Hilbert spaceD(A(q) ) endowed with the graph Hilbert space norm associated with the operatorA(q) . It follows (Curtain and Salamon, [12]) that Z =D(A(q) ),!H which means that Z is densely and continuously embedded in H. Dene Z to be the dual space of Z , we obtain that H ,!Z SinceA(q) is a bounded operator from Z toH, by duality, A(q) extends to a bounded operator ~ A(q) :H!Z. To nish this section, we enumerate all three operators and the spaces they 26 act on so as to be able to distinguish them below. A(q) :V !V A(q) :D(A(q))!H ~ A(q) :H!Z 4.2 TheAbstractCauchyProblemwithUnbounded Input In this section, we will summarize what has been done on linear Cauchy problems with unbounded input. Then we can apply these results to obtain well-posedness for abstract linear parabolic systems with input on the boundary of the spatial domain in the next section. We consider the abstract Cauchy problem _ '(t) =A(q)'(t) +Bu(t) '(0) =' 0 2H y(t) =C'(t) (4.2.1) on Hilbert space H where A(q) :D(A(q))!H is the innitesimal generator of a C 0 semigroup S(t) on H and C2L (H;Y ). Firstly, we want to be precise about what we mean by an unbounded input operator. Since the operatorA(q) is densely dened, closed and the innitesimal generator of aC 0 semigroup onH, it has an adjoint operatorA(q) :D(A(q) ) H! H, which is densely dened and closed. Dene the space Z to be the Hilbert spaceD(A(q) ) endowed with the graph Hilbert space norm associated 27 with the operator A(q) . It follows (Curtain and Salamon, [12]) that Z =D(A(q) ),!H which means that Z is densely and continuously embedded in H. Dene Z to be the dual space of Z , we obtain that H ,!Z Since A(q) is a bounded operator from Z to H, by duality, A(q) extends to a bounded operator ~ A(q) :H!Z. The input operator B then has range in the larger spaceZ. Then we want to dene a strong solution to Cauchy problem (4.2.1) in the space Z. Denition 4.2.1. (Strong Solution) '(t)2C 1 ([0;T ];Z)\C([0;T ];H) is said to be a strong solution of (4.2.1) if '(t) satises (4.2.1). This means that '(t) ='(0) + Z t 0 [ ~ A(q)'(s) +Bu(s)]ds; 0tT However, we want the solution to be in the space H where the output oper- ator is dened. Therefore we make the following hypothesis. Assumption. B1 For every T > 0, there exists a constant b T > 0 such that R T 0 S(Ts)Bu(s)ds2H and j Z T 0 S(Ts)Bu(s)dsjb T jju(t)jj Lp([0;T];U) for every u(t)2L p ([0;T ];U) where 1p<1. Theorem 4.2.1. If the above hypothesis is satised and let ' 0 2H and u(t)2 28 L p ([0;T ];U) be given, then '(t) =S(t)' 0 + Z t 0 S(ts)Bu(s)ds2H (4.2.2) is the unique strong solution of (4.2.1). Proof. Assume there are two solutions of (4.2.1) corresponding to the same input u() and the same initial state ' 0 ,denoted by '(t) and (t). So the dierence (t) ='(t)(t) satises (t) = Z t 0 ~ A (s)ds; 0tT Since ~ A(q) : H ! Z is a bounded operator, we obtain that (t)2 Z is continuously dierentiable on [0;T ] and satises _ (t) = ~ A(q) (t), (0) = 0. From semigroup theory, we know that (t) = 0 for 0tT . Now let '(t) be given by (4.2.2). Moreover, let t;s 2 [0;T ] and dene u(t) := 0 for t = 2 [0;T ]. Then, we have j'(t)'(s)j jS(t)' 0 S(s)' 0 j +j Z t s S()B[u(t)u(s)]dj jS(t)' 0 S(s)' 0 j +b T [ Z t s jju(t)u(s)jj p U d] 1=p and hence '(t) is continuous. Moreover, it is well known that for every contin- uously dierentiable input u(t), the function '(t) =S(t)' 0 + Z t 0 S(ts)Bu(s)ds; 0tT is continuously dierentiable and satises (4.2.1). So (4.2.2) is satised in this case. In general, (4.2.2) is satised from the fact that both sides of this equation 29 depend continuously on u(t)2L p ([0;T ];U). 4.3 AbstractLinearInitialBoundaryValueProb- lems with Input and Output on the Bound- ary of the Spatial Domain In this section, we will present similar results from (Curtain and Salamon, [12]) in dierent spaces. Assume V ,! H ,! V is a Gelfand triple, and W ,!V with dense and continuous embedding. Let U =R and Y =R . Let Q be a compact metric space, assume that for each q2 Q, we have bounded operators (q)2L (W;H), (q)2L (W;U) which is assumed surjective, and C(q)2L (V;Y ). Further assume that N((q)) is dense in H for each q2 Q. Consider the following initial-boundary value system with unbounded input on the boundary, _ '(t) = (q)'(t);t> 0 (q)'(t) =u(t) '(0) =' 0 y(t) =C(q)'(t) (4.3.1) Denition 4.3.1. (Strong Solution) Given u(t)2 C([0;T ];U) and ' 0 2 W such that (q)' 0 = u(0). '(t)2 C([0;T ];W )\C 1 ([0;T ];H) satisfying (4.3.1) for every t2 [0;T ] is said to be a strong(or classical) solution of (4.3.1). From Section 4.1, we know that under Assumptions (A1)-(A3), a(;;q) de- nes a bounded linear operator A(q) : V ! V such thata('; ;q) =< A(q)'; > V ;V , where '; 2 V . Furthermore, it can be shown (Banks and Ito, [6], Tanabe, [38]) that A(q) restricted on the set D(A(q)) =f'2 V : A(q)'2Hg so that A(q) =A(q) onD(A(q)) is closed and densely dened on 30 H and is the innitesimal generator of an analytic semigroup S(t;q) on H. It can further be shown thatA(q) : V V ! V is the innitesimal generator of an analytic semigroup onV that is an extension of the semigroupS(t;q) on H(see Tanabe, [[38]]). We assume that (q) is surjective and that its null space, N((q)) =f'2 W : (q)' = 0g is dense in H. Then assume further that the operator A(q) : Dom(A(q))2H!H satises: (1) Dom(A(q)) =f'2V :A'2Hg =f'2W : (q)' = 0g (2) A(q) is the restriction of (q) to Dom(A(q)). It follows that for each q2Q, (A(q);D(A(q))) = ((q);N((q)): (4.3.2) Note that the above assumption is satised for the alcohol model we devel- oped earlier. We will give the proof in Chapter 9. In the following discussion, on occasion we suppress showing the explicit dependence on t, that is when we write ' it denotes '(t) at a given t. We can now construct the input operator B(q)2L (U;V ). Given u2 U, we can choose w2 W such that (q)w = u since (q) is surjective. We then dene B(q)u := (q)wA(q)w2V : B(q) is well dened because (q)w = 0 if and only ifA(q)w = (q)w. It is easy to obtain the following properties of B(q): (i) B(q)(q)' = (q)'A(q)', '2W , (ii) B(q) is a linear map, (iii)jjB(q)ujj V jj(q)A(q)jj L (W;V ) jjwjj W C 1 jj(q)A(q)jj L (W;V ) jjujj U and therefore B(q) :U!V is bounded. 31 Indeed, for each q2 Q, let + (q)2L (U;W ) be any right inverse of the surjection (q), then B(q)2L (U;V ) by B(q) = ((q)A(q)) + (q) where we regard A(q)2 (V;V ). Note that B(q) is linear and well-dened regardless of the choice of + (q). Indeed, for two choices + 1 (q) and + 2 (q), and any q2Q, Range ( + 1 (q) + 2 (q)) N((q)). SinceA(q) = (q) on N((q)), B 1 (q) =B 2 (q). Because (q) was assumed to be surjective, there exists a constant K 0 > 0 such that for every u2U, there exists x2W such that (q)x =u;jjxjj W K 0 jjujj U : BecauseA(q)2L (V;V ), we also haveA(q)2L (W;V ) by continuity of the embedding W ,!V , since then for some a;a 0 , jjA(q)xjj V ajjxjja 0 jjxjj W : Therefore if u2U, then there exists x2W such that jjB(q)ujj V jj(q)A(q)jj L (W;V ) jjxjj W K 0 jj(q)A(q)jj L (W;V ) jjujj U ; so B(q)2L (U;V ). Assume that ' is a strong solution of (4.3.1). Then because B(q)(q) = (q)A(q), we have that _ '(t) = (B(q)(q) +A(q))'(t) =A(q)'(t) +B(q)u(t): 32 Therefore if ' is a strong solution of (4.3.1), then it is a solution to the nonhomogeneous Cauchy problem with unbounded input _ '(t) =A(q)'(t) +B(q)u(t);t> 0 '(0) =' 0 y(t) =C(q)'(t) (4.3.3) Note that replacing (A(q);D(A(q))) with A(q)2 L (V;V ), we can also regard this as a nonhomogeneous Cauchy problem inV with unbounded output since we have assumed thatC(q) has domainV ; in this case'(t)2V in general. Denition 4.3.2. (Weak Solution) LetA(q)2L (V;V ) andB(q)2L (U;V ), and let ' 0 2 V and u2 L 2 ([0;T ];U) be given. Then '(t)2 C([0;T ];V )\ W 1;2 ([0;T ];V ) is a weak solution of (4.3.1) if it satises (4.3.3) for almost every t2 [0;T ]. Denition 4.3.3. Suppose we are given a nonhomogeneous Cauchy problem (4.3.3), and assume thatA(q) is the innitesimal generator of an analytic semi- group S(t;q) on V , H, and V . B(q) will be said to be strongly p-admissible forA(q) if for some t 1 > 0, p2 [1;1], q2Q, the range of the linear map B(t;q) :L p ([0;t 1 ];U)!V u! Z t1 0 S(t 1 s;q)B(q)u(s)ds is actually in V . Note that to check p-admissibility it is enough to check for u in a dense subset of L p ([0;t 1 ];U). Lemma 4.3.1. Let '2H and u2U satisfyA(q)' +B(q)u2H, then '2W; (q)' =u; (q)' =A(q)' +B(q)u: 33 Moreover, there exists a constant C such that jj'jj W C[j'j +jjujj U +jA(q)' +B(q)uj] for all '2H, u2U withA(q)' +B(q)u2H. Proof. For '2 H and u2 U satisfying A(q)' +B(q)u2 H, choose w2 W such that (q)w =u (recall (q) was assumed to be surjective). Then A(q)('w) =A(q)'+B(q)u(A(q)+B(q))u =A(q)'+B(q)u(q)w2H this implies that '2W and jj'jj W jjwjj W +jj'wjj W jjwjj W +C 1 [j'wj +j(A(q)'w)j] jjwjj W +C 1 [j'j +jwj +jj(q)jj L (W;H) jjwjj W +jA(q)' +B(q)uj] [1 +C 1 jjIjj L (W;H) +C 1 jj(q)jj L (W;H) ]jjwjj W +C 1 j'j +C 1 jA(q)' +B(q)uj C 0 [1 +C 1 jjIjj L (W;H) +C 1 jj(q)jj L (W;H) ]jjwjj W +C 1 j'j +C 1 jA(q)' +B(q)uj C[j'j +jjujj U +jA(q)' +B(q)uj] Since (q)' = (q)w =g, we can get (q)' =A(q)' +B(q)(q)' =A(q)' + B(q)u. Theorem 4.3.2. Suppose that operator A(q) is the innitesimal generator of a strongly continuous semigroup S(t;q) on space H and that (q) is onto. If B(q) and A(q) satisfy Assumption B1 in Section 4.2, for every ' 0 2 V , u2L 2 ([0;T ];U), there exists a weak solution'2C([0;T ];V )\W 1;2 ([0;T ];V ) 34 of (4.3.1) given by '(t) =S(t;q)' 0 + Z t 0 S(ts;q)B(q)u(s)ds2H; t2 [0;T ]; (4.3.4) Furthermore, the strong solution of (4.3.1) is also given by (4.3.4). Proof. Suppose '(t) is given by (4.3.4) with ' 0 2 W and u(t)2 C 1 ([0;T ];U) satisfying (q)' 0 = u(0). Then '(t)2 C([0;T ];V )\C 1 ([0;T ];V ) and for t2 [0;T ], _ '(t) =A(q)'(t) +B(q)u(t) =S(t;q)[A(q)' 0 +B(q)u(0)] + Z t 0 S(ts;q)B(q) _ u(s)ds =S(t;q)(q)' 0 + Z t 0 S(ts;q)B(q) _ u(s)ds Using strong 2-admissibility, this implies _ '(t)2C([0;T ];H), sup t2[0;T] j _ '(t)jjS(t;q)jjj(q)jj L (W;H) jj' 0 jj W +M B (q)jj _ ujj L2([0;T];U) sup t2[0;T] M(q)e !(q)t jj(q)jj L (W;H) jj' 0 jj W +cM B (q)jj _ ujj L2([0;T];U) where we have usedjjcjjjj. Applying Lemma 4.3.1 to the termA(q)'(t) + B(q)u(t) = _ '(t)2H, we see that'(t)2C([0;T ];W ). Therefore'(t) is a strong solution of (4.3.1). BecauseC 1 ([0;T ];U) is densely and continuously embedded in C([0;T ];U), the result follows for arbitrary u(t)2C([0;T ];U). 35 Chapter 5 Abstract Semilinear Initial-Boundary Value Problems with Input on the Boundary of the Domain We want to extend the results from the last section to an abstract Initial- Boundary Value Problem with a nonlinear perturbation, _ '(t) = (q)'(t) +F (t;'(t);u(t);q);t> 0 (q)'(t) =G(t;'(t);u(t);q) '(0) =' 0 (q) y(t) =C(q)'(t) (5.0.1) where spaces H;V;W;U;Y;Q and operators (q); (q);C(q) are dened the same as in Section 4.3. Furthermore, we assume F (;;;q) : RWU! H and G(;;;q) :RWU!U. Firstly, we want to state what we mean by a strong solution to the above system precisely. Denition 5.0.1. (Strong Solution) Let u2 C[0;T ;U] and ' 0 2 W satisfy (q)' 0 =G(0;' 0 ;u(0)). Then a function'(t)2C[0;T ;W ]\C 1 [0;T ;H] is said to be a strong solution of (5.0.1) if (5.0.1) is satised for every t2 [0;T ]. In order to transform this abstract boundary value problem into a Cauchy problem, we make the same assumption as in the last section, (A(q);D(A(q))) = ((q);N((q))): 36 Therefore, we can also construct B(q) the same way as in Section 4.3, B(q) = ((q)A(q)) + (q): To be more precise, we will dene weak solution as below and state the equivalence theorem for strong solution and weak solutions. Denition 5.0.2. (Weak Solution) Let the operators A(q)2 L (V;V ) and B(q)2L (U;V ) be dened as above,' 0 2H be given, then'(t)2C(0;T ;H)\ W 1;2 [0;T ;V ] is said to be a weak solution of system(5.0.1) if _ '(t) =A(q)'(t) +F (t;'(t);u(t)) +B(q)G(t;'(t);u(t)) '(0) =' 0 (5.0.2) is satised in V for almost every t2 [0;T ]. It is a well-known semigroup theoretic result that the solution of (5.0.2) is given by '(t) =S(t;q)' 0 + Z t 0 S(ts;q)[F (s;'(s);u(s)) +BG(s;'(s);u(s))]ds (5.0.3) where S(t;q) is the semigroup generated by A(q) on H. In order to make sure that the solution is in the space H, we make the following hypothesis similar to the hypothesis in Section 4.2. Assumption. B2 For every T > 0, there exists a constant b T > 0 such that R T 0 S(Ts)BG(s;'(s);u(s))ds2H and j Z T 0 S(Ts)BG(s;'(s);u(s))dsjb T jju(t)jj Lp([0;T];U) for every u(t)2L p ([0;T ];U) where 1p<1. 37 Under the above assumption, we obtain that the weak solution of (5.0.2) '(t) =S(t;q)' 0 + Z t 0 S(ts;q)[F (s;'(s);u(s)) +BG(s;'(s);u(s))]ds2H In the discussion to follow, we will suppress showing all the arguments of the functions F (t;'(t);u(t)), G(t;'(t);u(t)) and refer to them as simply F;G respectively when no confusion will result. Theorem 5.0.1. A strong solutions and a weak solutions to System(5.0.1) are equivalent in the sense that a strong solution to system(5.0.1) is also a weak solution to system(5.0.1) as in denition 5.0.2, and every weak solution to system(5.0.1) which satises the regularity conditions in Denition(5) is also a strong solution to system(5.0.1). Proof. Suppose '2C([0;T ];W )\C 1 ([0;T ];H) is a strong solution, that is, _ '(t) = (q)'(t) +F (q)'(t) =G '(0) =' 0 We can then get _ '(t) = (q)'(t) +F +B(q)(G (q)'(t)) =A(q)'(t) +F +B(q)G and '(0) =' 0 which is (5.0.2). Conversely, suppose '(t)2 C([0;T ];W )\C 1 ([0;T ];H) is a weak solution to system(5.0.1) with ' 0 2W satisfying (q)' 0 =G(0;' 0 ;u(0)). Since '(t)2 C 1 ([0;T ];H), it follows thatA(q)'(t) +B(q)G = _ '(t)F 2 H. By Lemma 4.3.1, we have that '(t) satises (q)'(t) = G and (q)'(t) = A(q)'(t) + 38 B(q)G, which means that '(t) is also a strong solution to system(5.0.1). What's more, if we make the following assumptions onF ,G,C(q) and' 0 , we will be able to prove the existence and uniqueness of a weak solution(5.0.3) as in Denition 5. In the following discussion,jj;jjjj denote norms onH andV respectively. Assumption. C1 The mapping (t;';u)!G(t;';u) is continuous on RH U. Assumption. C2 There is a continuous mappingt!M(t) such thatG(t; 0;u) = M(t)u for all (t;u)2 RU. Furthermore, there exist m i 2L loc 1 ;i = 1; 2 such that jjG(t;';u;q)jj U fm 1 (t) +m 2 (t)jugj'j +jM(t)jjjujj U for all t;u and '. Assumption. C3 There is a continuous function g :RH!R, jjG(t;';u;q)G(t;';v;q)jj U g(t;')jjuvjj U for all (t;')2RH and u;v2U. Assumption. C4 For any > 0;q; ~ q2 Q, there exists q > 0 such that if jjq ~ qjj R p < q , then we have jjG(t;'(t;q);u(t);q)G(t;'(t; ~ q);u(t); ~ q)jj U <: Similarly, we make the following assumptions on F (;;) :RHU!H: Assumption. D1 The mapping (t;';u)!F (t;';u) is continuous on RH U. 39 Assumption. D2 There is a continuous mappingt!N(t) such thatF (t; 0;u) = N(t)u for all (t;u)2 RU. Furthermore, there exist n i 2 L loc 1 ;i = 1; 2 such that jF (t;';u;q)jfn 1 (t) +n 2 (t)jugj'j +jN(t)jjjujj U for all t;u and '. Assumption. D3 There is a continuous function f :RH!R jF (t;';u;q)F (t;';v;q)jf(t;')jjuvjj U for all (t;')2RH and u;v2U. Assumption. D4 For any > 0;q; ~ q2 Q, there exists 0 q > 0 such that if jjq ~ qjj R p < 0 q , then we have jF (t;'(t;q);u(t);q)F (t;'(t; ~ q);u(t); ~ q)j<: In addition, we make the following assumptions on ' 0 and C(q): Assumption. E For any and q; ~ q2Q, there exists ^ q > 0, ifjjq ~ qjj< ^ q , then we have j' 0 (q)' 0 (~ q)j<: Assumption. F For any and q; ~ q2 Q, there exists ~ q > 0, ifjjq ~ qjj < ~ q , then we have jjC(q)'C(~ q)'jj Y <: Theorem 5.0.2. UnderAssumptions (B2), (C1)-(C4), (D1)-(D4), (E), and (F), Equation(5.0.3) denes, for each' 0 2H, a unique functiont!'(t). Proof. Existence: 40 Suppose A(q) is the innitesimal generator of a strongly continuous semi- group S(t) on space H. Dene the Picard iteratesf' k (t)g by ' 0 (t) =S(t)' 0 ' k (t) =S(t)' 0 + Z t 0 S(ts)[F (s;' k1 (s);u(s)) +BG(s;' k1 (s);u(s))]ds (5.0.4) for k = 0; 1; 2;:::, and for t2 [0;t 1 ]. Clearly, the iterates are well dened, and an inductive argument employing the strong continuity offS(t)g and Assumption (B2) yields easily that ' k 2 C([0;t 1 ];H) for k = 0; 1; 2 . Using Assumptions (C2) and (D2) in the second equation of (5.0.4) j' k (t)jMe t fj' 0 j + Z t 0 jBjfm 1 (s) +m 2 (s)ju(s)gj' k1 (s)j +jM(s)jju(s)j +fn 1 (s) +n 2 (s)ju(s)gj' k1 (s)j +jN(s)jju(s)jdsg Me t fj' 0 j + + Z t 0 ~ U (s)j' k1 (s)jdsg where and ~ U (t)2L 2 ([0;T ]) is independent of k. A simple inductive argument using this inequality can be made from which one obtains j k (t)jp 0 Me t1 exp[H(t 1 )] where p 0 = sup t2[0;t1] fjS(t)' 0 j + g H(t) =Me t1 Z t 0 ~ U ()d Thus, it follows thatf' k (t)g is bounded in C([0;t 1 ];H). 41 From our assumptions, we have that jF (s;' k (s);u(s))F (s;' k1 (s);u(s))jfm 1 (s) +m 2 (s)ju(s)gj' k (s)' k1 (s)j for all k = 1; 2; . This can then be used to show that ' k (t) is in fact a Cauchy sequence in C([0;T ];H). Using completely straightforward arguments, one obtains j' k (t)' k1 (t)j C fMe t g k+1 c 0 [(t 1 ) k =k!] where c 0 = Z T 0 jF (s;' 0 (s);u(s))jds (t) = Z t 0 fm 1 (s) +m 2 (s)ju(s)gds Then, the usual arguments yield that the limit of the Cauchy sequencef' k g is the desired solution to system (5.0.1). Uniqueness: For any two continuous solutions ' and to system (5.0.1) on [0;T ], using Assumptions (C2) and (D2), for '(t);(t) lie in a bounded subsetD, j'(t)(t)j Z t 0 Me (ts) [fn 1 (s) +n 2 (s)ju(s)gj'(s)(s)j +Bfm 1 (s) +m 2 (s)ju(s)gj'(s)(s)j]ds Z t 0 U (s)j'(s)(s)jds whereU 2L 1 ([0;t]). By Gronwall's Inequality, we nd that j'(t)(t)jj'(0)(0)j exp( Z t 0 U (s)ds) = 0; 42 and consequently that '(t) =(t). 43 Chapter 6 Sampled Time Formulation For simplicity, when no confusion will result, in our discussions in this sec- tion and the next, we will suppress showing the explicit dependence on the parameters q in the operators and nonlinear perturbations. In the following discussion,jj denotes norm on H. Let the sampling time > 0 be given and u i be a zero-order hold input u i =u(i), i = 0; 1; 2; . Based on the solution to the abstract parabolic system given in (5.0.3), dene sequencef' i g: ' 0 ='(0) ' i+1 =S(;q)' i + Z 0 S(;q)dF (i;' i ;u i ) + Z 0 S(;q)B(q)dG(i;' i ;u i ) y i =C(q)' i (6.0.1) SinceS(t;q) is an analytic semigroup onV;H andV , it follows thatS(;q)2 L (V;V ) and R 0 S(;q)B(q)d2L (U;V ). So for simplicity, we can write the above discrete innite dimensional system as ' 0 ='(0) ' i+1 =S(;q)' i + ^ A(q)F (i;' i ;u i ) + ^ B(q)G(i;' i ;u i ) y i =C(q)' i where ^ A(q) = R 0 S(;q)d and ^ B(q) = R 0 S(;q)B(q)d. Theorem 6.0.1. As ! 0,j'(i)' i j! 0 for all i = 0; 1; 2; . 44 Proof. j'((i + 1))' i+1 jjS(;q)jj'(i)' i j +j Z 0 S(;q)(F (i +;'(i +);u(i +))F (i;' i ;u i ))dj +j Z 0 S(;q)(B(q)G(i +;'(i +);u(i +))B(q)G(i;' i ;u i ))dj :=I 1 +I 2 +I 3 We can rst look at the second term, I 2 Z 0 jS(;q)jjF (i +;'(i +);u(i +))F (i;'(i +);u(i +))jd + Z 0 jS(;q)jjF (i;'(i +);u(i +))F (i;'(i);u(i +))jd + Z 0 jS(;q)jjF ('(i);u(i +))F (i;' i ;u(i +))jd + Z 0 jS(;q)jjF (' i ;u(i +))F (i;' i ;u i )jd By Assumptions (C1) to (C3), we obtain I 2 Z 0 jS(;q)jjF (i +;'(i +);u(i +))F (i;'(i +);u(i +))jd + Z 0 Me fm 1 () +m 2 ()ju(i +)gj'(i +)'(i)jd + Z 0 Me fm 0 1 () +m 0 2 ()ju(i +)gj'(i)' i jd + Z 0 Me jg(' i )jju(i +)u i jd U ()j'(i)' i j + Z 0 a()j'(i +)'(i)jd + Z 0 Me jg(' i )jju(i +)u i jd where U () = Me + R 0 Me fm 1 () + m 2 ()ju(i + )gd and a() = Me fm 0 1 () +m 0 2 ()ju(i +)g. 45 By the strong continuity ofS(t) and continuity ofF (;;), we can infer that '(t)2 C([0;T ];V ). Therefore, the second term and the third term go to 0 as ! 0. Also, the third term goes to 0 as ! 0 since u(t)2C([0;T ];R n ). SinceU () is bounded, we can prove iteratively thatj'((i+1))' i+1 j! 0. Similarly, I 3 ! 0 Under Assumption (A2), with possinly a change of variable, we may assume without loss of generality, that the operatorA is invertible with bounded inverse. In this case, it is possible to eliminate the integrals in the above expression, equation(6.0.1) and obtain that ' i+1 =S(;q)' i + Z 0 S(;q)dF (i;' i ;u i ) + Z 0 S(;q)dBG(i;' i ;u i ) =S(;q)' i + (S(;q)I)A 1 F (i;' i ;u i ) + (S(;q)I)A 1 (A) + G(i;' i ;u i ) =S(;q)' i + Z 0 S(;q)dF (i;' i ;u i ) + Z 0 S(;q)d(q) + (q)G(i;' i ;u i ) (S(;q)I) + (q)G(i;' i ;u i ) So with the assumption that ' 0 2 V and F (i;' i ;u i )2 V , we have that ' i 2V , for i = 1; 2; 3; . 46 Chapter 7 Finite Dimensional Approximation and Convergence For each N = 1; 2; , let H N be a nite dimensional subspace of V with basisf N j g N j=1 with the property that for every 2V , there exists N 2H N such that lim N!1 N = in V . Let P N : H! H N denote the orthogonal projection of H onto H N . Dene A N (q) : H N ! H N to be a bounded linear operator onH N and letS N (t;q) be aC 0 semigroup generated byA N (q) for all t> 0, ' i+1 =S(;q)' i + Z 0 S(;q)dF (i;' i ;u i ) + Z 0 S(;q)d(q) + (q)G(i;' i ;u i ) (S(;q)I) + (q)G(i;' i ;u i ) (7.0.1) ' N i+1 =S N (;q)P N ' N i + Z 0 S N (;q)dP N F (i;' N i ;u i ) + Z 0 S N (;q)dP N (q) + (q)G(i;' N i ;u i ) (S N (;q)I)P N + (q)G(i;' N i ;u i ) (7.0.2) or ' i+1 =S(;q)' i + ^ A(q)F (i;' i ;u i ) + ^ B(q)G(i;' i ;u i ) ' N i+1 =S N (;q)P N ' N i + ^ A N (q)F (i;' N i ;u i ) + ^ B N (q)G(i;' N i ;u i ) 47 where ^ A N (q) = R 0 S N (;q)dP N 2L (V;H N ) and ^ B N (q) = R 0 S N (;q)dP N (q) + (q) (S N (;q)I)P N + (q)2L (U;H N ). Lemma 7.0.1. Suppose thatf i g 1 i=0 andf i g 1 i=0 are sequences of nonnegative real numbers and thatf i g 1 i=0 is a sequence of real numbers which satisfy n n + n1 k=0 k k ; n = 1; 2;::: Then we have that n n + n1 j=0 j j (exp( n1 k=j+1 k )); n = 1; 2;::: If in addition j = 0, j = 0; 1; 2;:::, then n exp( n1 k=0 k ) n = 1; 2;:::: Proof. The proof, which is not dicult, has been omitted but can be found in [26]. Theorem 7.0.2. Given ' i and ' N i are given by equations (7.0.1) and (7.0.2). Suppose (1) For every sequencefq N g 1 N=1 2 Q with lim N!1 q N = q, there exists 2(A(q))\(A N (q N )) with Re>! such that lim N!1 R (A N (q N ))P N = R (A(q)) in V , for each 2 V , where (A) denotes the resolvent set of a linear operator A, R (A) = (A) 1 denotes the resolvent operator of a linear operator A at 2(A). (2) There exist constants M;, such thatjS(t)j Me t ,jS N (t)j Me t , for N = 1; 2;::: Then as N!1,jj' N i ' i jj V , for i = 0; 1; 2;::: Proof. Let K 0 denote the uniform bound on the operators S N (t) for all N 48 suciently large, K 1 denote the uniform bound on the operators R 0 S N ()d, jj + Gjj V C 1 jjGjj U andjj + Gjj W C 2 jjGjj U for all N suciently large. jj' i ' N i jj V jjS(i)' 0 S N (i)P N ' 0 jj V + i1 k=0 jjS N (k)jj V jj Z 0 S N ()djj V jjP N F (k;' k ;u k )P N F (k;' N k ;u k )jj V + i1 k=0 jjS(k) Z 0 S()dF (k;' k ;u k )S N (k) Z 0 S N ()dP N F (k;' k ;u k )jj V + i1 k=0 jjS N (k)jjjj Z 0 S N ()djjjjP N G(k;' k ;u k )P N + G(k;' N k ;u k )jj V + i1 k=0 jjS(k) Z 0 S()d + G(k;' k ;u k )S N (k) Z 0 S N ()dP N G(k;' k ;u k )jj V + i1 k=0 jjS N (k)jj(S N ()I)jjP N + G(k;' k ;u k )P N + G(k;' N k ;u k )jj V + i1 k=0 jjS(k)(S N ()I) + G(k;' k ;u k )S N (k)(S N ()I)P N + G(k;' k ;u k )jj V + i1 k=0 K 0 K 1 [(C 1 +C 2 )fm 1 +m 2 jjug +fn 1 +n 2 jjug]jj' k ' N k jj V + for allN suciently large. If we now apply Lemma 7.0.1 to the above inequality, jj' i ' N i jj V 2 exp( ) for all N suciently large, where = i1 k=0 K 0 K 1 [(C 1 +C 2 )fm 1 +m 2 jug +fn 1 +n 2 jug] and the theorem is proven. 49 Chapter 8 Parameter Estimation In the case of parameter estimation, we assume that training data isf^ y i ; ^ u i g L i=1 , where ^ y i , ^ u i are sampled from continuous time signals ^ y and ^ u. Moreover, to implement numerical analysis, we will add a regularization termjjqq 0 jj 2 with weight ! that forces optimal parameters to be close to a given set of values q 0 . Find q 2Q which minimizes the quadratic J(q) = L i=1 jjy(i;q) ^ y i jj 2 +!jjqq 0 jj 2 where y(i;q) is given by equation (3.0.1). Find q N 2Q which minimizes the discrete quadratic form J N (q) = L i=1 jjy N i (q) ^ y i jj 2 +!jjqq 0 jj 2 where y N i (q) =C(q)' N i (q) and ' N i (q) is given by (7.0.2). Theorem 8.0.1. Given ' i (q) and ' N i (q N ) are given by equations (7.0.1) and (7.0.2). Suppose (1) For every sequencefq N g 1 N=1 2 Q with lim N!1 q N = q, there exists 2(A(q))\(A N (q N )) with Re>! such that lim N!1 R (A N (q N ))P N = R (A(q)) in V , for each 2 V , where (A) denotes the resolvent set of a linear operator A, R (A) = (A) 1 denotes the resolvent operator of a linear operator A at 2(A). (2) There exist constants M;, such thatjS(t)j Me t ,jS N (t)j Me t , for N = 1; 2;::: 50 Then as N!1,jj' N i (q N )' i (q)jj V , for i = 0; 1; 2;::: Proof. Denote R 0 S(;q)d as B(q) and R 0 S N (;q)d as B N (q), jj' N i (q N )' i (q)jj V jj' N i (q N )' N i (q)jj V +jj' N i (q)' i (q)jj V jj' N i (q)' i (q)jj V +jjS N (i;q N )P N ' 0 (q N )S N (i;q)P N ' 0 (q)jj V + i1 k=0 jjS N (k;q N )B N (q N )P N F (k;' N k (q N );u k ;q N ) S N (k;q)B N (q)P N F (k;' N k (q);u k ;q)jj V I 1 +I 2 +I 3 According to Theorem 7.0.2 , I 1 for all N suciently large. I 2 jjS N (i;q N )P N ' 0 (q N )S N (i;q N )P N ' 0 (q)jj V +jjS N (i;q N )P N ' 0 (q)S N (i;q)P N ' 0 (q)jj V Me t jj' 0 (q N )' 0 (q)jj V +jjS N (i;q N )P N S N (i;q)P N jjjj' 0 (q)jj V By the strong continuity of S(t) and continuity of ' 0 , I 2 ! 0 I 3 i1 k=0 jjS N (k;q N )B N (q N )P N F (k;' N k (q N );u k ;q N ) S N (k;q)B N (q)P N F (k;' N k (q);u k ;q)jj V i1 k=0 jjS N (k;q N )B N (q N )P N F (k;' N k (q N );u k ;q N ) S N (k;q N )B N (q N )P N F (k;' N k (q);u k ;q)jj V + i1 k=0 jjS N (k;q N )B N (q N )P N F (k;' N k (q);u k ;q) S N (k;q)B N (q)P N F (k;' N k (q);u k ;q)jj V By the strong continuity of S(t) and hypothesis (B4), I 3 ! 0 51 Due to the compactness ofQ there exists a subsequencef q Nj g 1 j=1 2f q N g 1 N=1 and q such that lim j!1 q Nj = q. We will prove that q solves the original parameter estimation problem. J( q) = L i=1 jy(i; q) ^ y i j 2 = L i=1 jC( q)' i ( q) ^ y i j 2 = L i=1 jC( lim j!1 q Nj ) lim j!1 ' i;Nj ( q Nj ) ^ y i j 2 = lim j!1 L i=1 jC( q Nj )' i;Nj ( q Nj ) ^ y i j 2 = lim j!1 J Nj ( q Nj ) lim j!1 J Nj (q) = L i=1 jC(q) lim j!1 ' i;Nj (q) ^ y i j 2 = L i=1 jC(q)' i (q) ^ y i j 2 =J(q) 52 Chapter 9 Application to our Alcohol Model In order to apply the theory developed in the previous chapters to our cou- pled ODE and PDE model for alcohol transport and metabolism we developed in Chapter 3, we must rst dene and identify the operatorsA(q),A(q), spaces H;V ;V;D(A(q)) and functions F ();G() for the given system below. @x @t (t;) =D @ 2 x @ 2 (t;); 0l;t> 0 dw(t) dt =A 0 w(t) Vw(t) K +w(t) +bu(t);t> 0 D @' @ (t; 0) ='(t; 0) D @' @ (t;l) =w(t) '(0;) = 0; 0l w(0) = 0 y(t) =r'(t; 0) where D;;r2R, A 0 2R nn , w()2C 1 ([0;T ];R n1 ) and 2R 1n . As the above system shows, we are solving a semi-linear parabolic system with homogeneous boundary condition, and we can identify functions F (t;x(t);w(t);u(t)) = Vw(t) K +w(t) +bu(t) and G(t;x(t);w(t);u(t)) = 0 . 53 We denote 2 6 4 x w 3 7 5 as ^ x, 2 6 4 ' r 3 7 5 as ^ in the following discussion. Let H := L 2 ([0;l])R, V := H 1 ([0;l])R and W := H 2 ([0;l])R, and their inner products as < ^ x; ^ > H = R l 0 x d +wr, < ^ x; ^ > V = R l 0 x d + R l 0 x 0 0 d +wr and < ^ x; ^ > W = R l 0 x d + R l 0 x 0 0 d + R l 0 x 00 00 d +wr, where ^ x = (x;w) and ^ = ( ;r). Denote H 1 0 ([0;l]) := f j 2 H 1 ([0;l]); (0) = 0; (l) = 0g, H 1 L ([0;l]) :=f j 2 H 1 ([0;l]); (0) = 0g and H 1 0 (0;l) :=f j 2 H 1 ([0;l]); (l) = 0g and inner product <;> without specifying space will be L 2 inner product <x; >= R l 0 x d. We can rst get the sesquilinear form a(^ x; ^ ;q) =D<' 0 ; 0 >w'(1) +'(0) (0)w T A T 0 r from Chapter 3. ThereforeA is dened by <A^ x; ^ > V ;V =D<x 0 ; 0 > +w (1) (0) (0) +w T A T 0 r where<;> V ;V now denotes the natural extension of the H innerproduct to the duality pairing between V and V and where we have suppressed showing explicit dependence on t wherever possible. As dened at the beginning of this section, we have D(A(q)) =f'2V :A(q)'2Hg and A(q) :D(A(q))H!H is dened by A(q) =A(q) onD(A(q)). On the other hand, we dene ^ D =f^ x2H 2 (0;l)R n :Dx 0 xj 0 = 0;Dx 0 wj l = 0g 54 and dene an operator ^ A : ^ DH!H as ^ A^ x = 2 6 4 D @ 2 @ 2 0 0 A 0 3 7 5 2 6 4 x w 3 7 5: (9.0.1) Theorem 9.0.1. A = ^ A Proof. In order to prove this theorem, we basically need to prove (1) ^ A A and (2) A ^ A. (1) To prove ^ AA, we need to prove that ^ DD(A(q)) and ^ A^ x =A^ x for ^ x2 ^ D. First we will show that ^ DD(A(q)), let ^ x2 ^ D and ^ 2V < ^ A^ x; ^ > H =<D @ 2 @ 2 x; > +w T A T 0 r =D @x @ j 1 0 <D @x @ ; 0 > +w T A T 0 r =w (1)x(0) (0)<D @x @ ; 0 > +w T A T 0 r =<A^ x; ^ > H =) if ^ x2 ^ D, thenA^ x = 2 6 4 @ 2 x @ 2 A 0 w 3 7 52L 2 ([0;l])R n =H, sincex2H 2 ([0;l]) =) if ^ x2 ^ D, then ^ x2D(A(q)) and for ^ x2 ^ D, ^ A^ x =A^ x =A^ x =) ^ AA. (2) Now to show that A ^ A, we will showD(A(q)) ^ D. Note if this is true, then ^ x2D(A(q))) ^ x2 ^ D and ^ A^ x =A^ x = A^ x, which would complete the proof. Let ^ x2D(A(q)), and letA^ x = ^ = 2 6 4 v 3 7 52H. 55 Then for ^ 2V , < ^ ; ^ > H =<A^ x; ^ > H =D Z l 0 x 0 0 d +w (l)x(0) (0) +w T A T 0 r <; > +v T r =D Z l 0 x 0 0 d +w (l)x(0) (0) +w T A T 0 r Z l d j l 0 Z l 0 Z l d 0 d +v T r =D Z l 0 x 0 0 d +w (l)x(0) (0) +w T A T 0 r Let 2H 1 0 ([0;l]) and r = 0, we can get Z 1 0 Z 1 d 0 d =D Z 1 0 x 0 0 d: Since H 1 0 ([0;l]) = L 2 ([0;l]), we can get R 1 d = Dx 0 in L 2 ([0;l]), which implies Dx 0 2H 1 ([0;l]), so that x2H 2 ([0;l]) and =D @ 2 x @ 2 . Now for ^ 2V , Z l 0 D @ 2 x @ 2 d +v T r = Z l 0 d +v T r =D Z l 0 x 0 0 d +w (l)x(0) (0) +w T A T 0 r )Dx 0 j l 0 D Z l 0 x 0 0 d +v T r =D Z l 0 x 0 0 d +w (l)x(0) (0) +w T A T 0 r )Dx 0 j l 0 +v T r =w (l)x(0) (0) +w T A T 0 r: Then let 2H 1 L ([0;l]) and r = 0, we can get Dx 0 (l) =w: Let 2H R ([0;l]) and r = 0, we can get Dx 0 (0) =x(0): 56 Now we can get v =A 0 w and we come to the conclusion that ^ x2 ^ D. Therefore, the abstract semi-linear parabolic system we are solving for the alcohol model is ' i+1 =S(;q)' i + ^ B(q)F (i;' i ;u i ;q) ' 0 = 0 y i =C(q)' i ;i = 1; 2; ;L where S(t;q) is a semigroup generated by ^ A(q) dened in (9.0.1) and ^ B(q) = R 0 S(;q)d. And the corresponding nite dimensional approximation system is ' N i+1 =S N (;q)' N i + ^ B N (q)F (i;' N i ;u i ;q) ' N 0 = 0 y N i =C(q)' N i ;i = 1; 2; ;L where S N (t;q) is a semigroup generated by ^ A N (q) dened in Chapter 7 and ^ B N (q) = R 0 S N (;q)d. 9.1 Estimating Blood or Breath Alcohol Con- centration from Biosensor Measured Trans- dermal Alcohol Concentration In this section, we will concentrate on applying the above theorems we de- veloped for the abstract system to the blood alcohol concentration estimation from biosensor. At rst, We need to check that alcohol concentration model in 57 Chapter 3 satises all the assumptions. (A1) Continuity: Proof. ja(; ;q)a(; ;e q)jjq 1 e q 1 j Z l 0 jx 0 z 0 jd +jx(0)z(0)j +jq 4 e q 4 jjwjjz(1)j +jw T (A T e A T )rj jq 1 e q 1 jjjx 0 jj L2 jjz 0 jj L2 +jx(0)jjz(0)j +jq 4 e q 4 jjwjjz(1)j +jwjjA T e A T j 1 jrj d(q;e q)jjjj V jj jj V (A2) Coercivity: Proof. a(;;q) + 0 jjjj 2 H = 0 jjjj H +q 1 Z l 0 (x 0 ) 2 d +x 2 (0)w T A T wq 4 wx(1) We now work on the last summand: q 4 wx(1)q 4 Cwjjxjj 2 by Sobolev's Inequality 58 Then for some large 0 , a(;;q) + 0 jjjj 2 H 0 Z l 0 x 2 d + 0 w 2 +q 1 Z l 0 (x 0 ) 2 d +x 2 (0)w T A T wq 4 Cwjjxjj 2 c 1 jjjj V where c 1 =min 0 ;q 1 ; 0 A T ;q 4 (A3) Boundedness: Proof. ja(; ;q)jjq 1 j Z l 0 jx 0 z 0 jd +jx(0)z(0)j +jq 4 jjwjjz(1)j +jw T A T rj jq 1 jjjx 0 jj L2 jjz 0 jj L2 d +jx(0)jjz(0)j +jq 4 jjwjjz(1)j +jwjjA T j 1 jrj c 2 jjjj V jj jj V where c 2 =maxfjq 1 j;jq 4 j;jA T j 1 ; 1g (D1) Proof. For ' = (x;w) and = ( ;r) j kw M +w +bu ( kr M +r +bu)jj kw(M +r)kr(M +w) (M +w)(M +r) j j kM (M +w)(M +r) jjwrj (D2) 59 Proof. j kw M +w +bujjkjjwj +jbjjuj (D3) Proof. For u and v j kw M +w +bu ( kw M +w +bv)jjbjjuvj Since G(t;x(t);w(t);u(t)) = 0, so Assumptions (C1)-(C4) are automati- cally true. Assumptions (B2), (D4), (E) and (F) are obvious to prove, since F () and ' 0 are continuous with respect to k and M. 9.2 The Adjoint Method In order to implement numerical experiments, we will use the fully discretized model. What's more, to implement numerical analysis, we will add a regular- ization termjjqq 0 jj 2 with weight! that forces optimal parameters to be close to a given set of values q 0 . Minimize over all q2Q the least squares functional J N :Q!R given by J N (q) = L X i=1 jjy N i (q) ^ y i jj 2 +!jjqq 0 jj 2 (9.2.1) 60 subject to the abstract discrete system below ' N i+1 =S N (;q)' N i + ^ B N (q)F (i;' N i ;u i ;q) ' N 0 = 0 y N i (q) =C N (q)' N i ;i = 1; 2; ;L where S N (t;q) is a semigroup generated by ^ A N (q) dened in Chapter 7 and ^ B N (q) = R 0 S N (;q)d. Now we need to solve the above discretized optimization problem. In order to use steepest descent or quasi-Newton method to search for local minima, we need to compute the gradientrJ N . Using nite dierences to compute the gradient would be computational expensive. Since we already have a fully discrete approximation, we apply the adjoint method to compute the gradient rJ N analytically without truncation error. We will derive adjoint method here, @J N (q) @q = 2 L i=1 [C(q)' N i (q) ^ y i ] T [C(q) @' N i (q) @q + @C(q) @q ] + 2!jqq 0 j @' N i+1 (q) @q =A N (q) @' N i (q) @q + @A N (q) @q ' N i (q) + @B N (' N i ;q) @' N i @' N i (q) @q + @B N (' N i ;q) @q where we denote B N (' N i ;q) = ^ B N (q)F (i;' N i ;u i ) We dene the co-states z i , and the adjoint equation is given by z i1 = [A N (q) + @B N (' N i1 ;q) @' N i1 ] T z i + i1 (q) z L = 2!jqq 0 j where i (q) = 2C T (q)e i (q) with e i (q) = [C(q)' N i (q) ^ y i ]. 61 It then follows that @J N (q) @q = L i=0 T i (q) @' N i (q) @q + 2 L i=0 e T i (q) @C(q) @q ' N i (q) = L i=0 [z i (A N (q) + @B N (' N i ;q) @' N i ) T z i+1 ] T @' N i (q) @q + 2 L i=0 e T i (q) @C(q) @q ' N i (q) = L i=0 z T i @' N i (q) @q L i=0 z T i+1 (A N (q) + @B N (' N i ;q) @' N i ) @' N i (q) @q + 2 L i=0 e T i (q) @C(q) @q ' N i (q) =z T 0 @' N 0 (q) @q + L i=1 z T i @' N i (q) @q L i=1 z T i (A N (q) + @B N (' N i1 ;q) @' N i1 ) @' N i1 (q) @q + 2 L i=0 e T i (q) @C(q) @q ' N i (q) =z T 0 @' N 0 (q) @q + L i=1 z T i (q)! i1 (q) + 2 L i=0 e T i (q) @C(q) @q ' N i (q) where ! i (q) = @A N (q) @q ' N i (q) + @B N (' N i ;q) @q . To summarize, to carry out adjoint method, we give Algorithm 1. Algorithm 1: Adjoint Method Input: S N (;q),B N (' N i ;q) Output:rJ N (q) 1 Update co-states Set z M (q) = 2!jqq 0 j. For i =L;:::; 1, z i1 = [A N (q) + @B N (' N i ;q) @' N i1 ] T z i1 + i1 (q) 2 Compute gradient rJ N (q) =z T 0 @' N 0 (q) @q + L i=1 z T i (q)! i1 (q) + 2 L i=0 e T i (q) @C(q) @q ' N i (q) where ! i (q) = @A N (q) @q ' N i (q) + @B N (' N i ;q) @q . 62 Chapter 10 Numerical Studies for the Deterministic Model 10.1 Deterministic Model We consider the model for alcohol metabolism and transport derived in the previous chapters given by @x @t (t;) =D @ 2 x @ 2 (t;); 0L;t> 0 dw(t) dt =aw(t) Vw(t) K +w(t) +bu(t);t> 0 D @x @ (t; 0) =x(t; 0) D @x @ (t;l) =w(t) x(0;) = 0; 0l w(0) = 0 y(t) =rx(t; 0): (10.1.1) As can be readily seen in (10.1.1), there are 7 parameters D;a;V;K;b; and r in the model that require estimation. Our objective in this section is to determine the optimal estimates of these parameters, which we will refer as the model parameters in the rest of this thesis to distinguish them from the drinking parameters which we will discuss later. In order to estimate the model parameters, we will use so called calibration input and corresponding output data sets,fu i ;y i g n i=1 . Not all of the parameters are independent and consequently they are not all identiable. As a result, we are able to simplify the above system into the 63 following system, @^ x @ ^ t ( ^ t; ^ ) = ^ D @ 2 ^ x @^ 2 ( ^ t; ^ ); 0 ^ 1; ^ t> 0 d ^ w( ^ t) d ^ t =a ^ w( ^ t) V ^ w( ^ t) K + ^ w( ^ t) +bu(t); ^ t> 0 ^ D @^ x @^ ( ^ t; 0) = ^ x( ^ t; 0) ^ D @^ x @^ ( ^ t; 1) = ^ w( ^ t) ^ x(0; ^ ) = 0; 0 ^ l ^ w(0) = 0 y(t) = ^ x( ^ t; 0) (10.1.2) where ^ x = x L , ^ t = t , ^ D = D L 2 , ^ a = a L , ^ = ar L , = L a and ^ x = rx( ^ t;L^ x), 0 ^ 1, ^ t 0. Consequently, in calibrating this revised model, rather than t the seven parameters D;a;V;K;b; and r in (10.1.1), it is sucient to t the ve parameters, ^ D;a;V;K and in (10.1.2). Therefore, the model we use in the parameter estimation problem is the following, @x @t (t;) =q 1 @ 2 x @ 2 (t;); 0 1;t> 0 dw(t) dt =q 2 w(t) q 3 w(t) q 4 +w(t) +u(t);t> 0 @x @ (t; 0) =x(t; 0) q 1 @x @ (t;l) =q 5 w(t) x(0;) = 0; 0l w(0) = 0 y(t) =x(t; 0) and the parameters are q = [q 1 ;q 2 ;q 3 ;q 4 ;q 5 ]. 64 10.2 Finite Dimensional Approximation In order to solve the optimization problem (9.2.1) numerically, we need to apply nite dimensional approximation to get operatorsA N (q) and ^ B N (q) that can now be represented as matrices. Based on the theory in Chapter 7, we discretize the spatial domain using linear B-splines and Galerkin approximation. ForN = 1; 2; , letf N j g N j=0 denote the set of standard linear B-splines on the interval [0; 1] dened in terms offj=Ng N j=0 , where N j is a "pup tent" function with height 1 and support of width 2=N, [(j 1)=N; (j + 1)=N]\ [0; 1]. Let f ^ N j g N j=0 =f( N j ; 1)g N j=0 . Let H := L 2 ([0;l])R, V := H 1 ([0;l])R and W := H 2 ([0;l])R, and their inner products as < ^ x; ^ > H = R l 0 x d +wr, < ^ x; ^ > V = R l 0 x d + R l 0 x 0 0 d +wr and < ^ x; ^ > W = R l 0 x d + R l 0 x 0 0 d + R l 0 x 00 00 d +wr, where ^ x = (x;w) and ^ = ( ;r). Denote H 1 0 ([0;l]) := f j 2 H 1 ([0;l]); (0) = 0; (l) = 0g, H 1 L ([0;l]) :=f j 2 H 1 ([0;l]); (0) = 0g and H 1 0 (0;l) :=f j 2 H 1 ([0;l]); (l) = 0g and inner product <;> without specifying space will be L 2 inner product <x; >= R l 0 x d. Let H N =spanf ^ N j g N j=0 and P N :H!H N denote the orthogonal projec- tion of H onto H N with respect to the H inner product. For N = 1; 2; , and q2 Q, A N (q) can be represented as a matrix given by, A N (q) = [< ^ N i ; ^ N j >] 1 [a( ^ N i ; ^ N j ;q)] where a( ^ N i ; ^ N j ;q) =q 1 < N 0 i ; N 0 j >w N j (1) + N i (0) N j (0)q 2 wr (10.2.1) To show how to compute operatorsA N (q) and ^ B N (q), we give the denition 65 of the following (N + 2) (N + 2) matrices, M = [< ^ N i ; ^ N j >] = 1 6N Tridiagf[1; 1;; 1; 0]; [2; 4; 4; ; 4; 2; 6N]; [1; 1;; 1; 0]g K = [< N 0 i ; N 0 j >] =NTridiag[1;1; ;1; 0]; [1; 2; 2; ; 2; 0]; [1;1; ;1; 0] L = [ N i (0) N j (0)] =diag[1; 0; 0; ; 0] R = [ N j (1)] =diag[0; 0; ; 0; 1] A =diag[0; 0; ; 0; 1]: From (10.2.1), we have A N (q) =M 1 (q 1 KL +Rq 2 A) and S N (;q) =e A N (q) . According to the denition of ^ B N (q) = R 0 S N (;q)d in Chapter 9, we have ^ B N (q) = (S N (;q)I)A N (q) 1 M 1 [0; 0; ; 0; 1] T Therefore, y N i (q) is given by ' N i+1 =S N (;q)' N i + ^ B N (q)F (i;' N i ;u i ;q) ' N 0 = 0 y N i (q) =C(q)' N i ;i = 1; 2; ;L 66 10.3 Experimental Data Description The experimental data we have is from one subject who wore a WrisTAS TM 7 alcohol biosensor for 18 days. The WrisTAS TM 7 was set to measure the alcohol concentration over the subject's wrist at 5-minute intervals. For each drinking episode, BrAC readings were taken every 30 minutes until the BrAC returned to 0.000. As shown in Figure 10.3.1, there are in total 11 individual drinking episodes recorded in this experiment. Figure 10.3.1: Experimental Data: BrAC and WrisTAS TM 7 alcohol biosensor TAC measurements. 67 10.4 Estimating the Model Parameters Based on One Drinking Episode In this section, we show how we solved the parameter estimation problem using codes written in MATLAB and present the numerical results. We assume that we have data from m drinking episodes and that the training data from the i-th drinking episode isf^ y ij ; ^ u ij g ni j=1 , i = 1; 2; ;m, where ^ y ij is the j-th TAC measurement, ^ u ij is the j-th measurement of the amount of alcohol ingestion. If we use the data from i-th drinking episode to estimate parameters q, i.e. as training data, then the optimization problem we need to solve is given by ^ q = arg min q2Q J N (q) = arg min q2Q m X i=1 ni X j=1 jy N i;j (q) ^ y ij j 2 (10.4.1) subject to the abstract discrete system ' N i;j+1 =S N (;q)' N i;j + ^ B N (q)F (j;' N i;j ; ^ u i;j ;q) ' N i;0 = 0 y N i;j (q) =C N (q)' N i;j ;j = 1; 2; ;n i ;i = 1; 2; ;m; (10.4.2) where S N (t;q) and ^ B N (q) are computed numerically as in Section 10.2. To solve the above optimization problem, we used the MATLAB routine FMINCON in the Optimization Toolbox. FMINCON requires two inputs: a function that evaluates the objective function and a function that evaluates the gradient of the objective function. Evaluation of the objective function is given by (10.4.2) and (10.4.1) denes the optimal set of parameters. The gradient 68 of the objective function can be computed using the adjoint method which is described in Section 9.2. Note that the optimization problem we solve here is unconstrained. And we ensured the positivity of the parameters by searching on the squares of the parameters rather than the parameters themselves. In order to show how well the model with the t parameters agrees with the data, we produced Table 10.4.1 with the optimal parameters and several metrics measuring the dierence between the tted BrAC (and TAC) and the experimentall observed BrAC (and TAC). Drinking Episode 1 2 3 4 5 6 q1 1.14E+00 1.95E+02 2.24E+02 5.59E-01 8.53E+00 4.13E+02 q2 2.77E-04 1.28E+00 7.74E-01 3.29E-02 2.25E-01 1.10E+00 q3 7.38E-01 8.29E-01 7.81E-01 8.05E-01 6.78E-01 1.22E+00 q4 1.60E-01 6.14E+03 1.80E+01 1.39E-01 5.98E+05 9.01E+03 q5 1.41E-01 4.26E+04 6.03E+01 2.01E-02 1.00E+06 2.61E+04 TAC L 2 3.42E-02 2.35E-02 6.15E-02 9.19E-02 1.15E-01 6.16E-02 TAC L 1 9.68E-03 6.08E-03 1.43E-02 1.22E-02 2.10E-02 1.05E-02 BrAC L 2 4.36E-02 2.67E-02 6.24E-02 9.07E-02 1.31E-01 4.04E-02 BrAC L 1 1.39E-02 8.82E-03 1.40E-02 1.86E-02 2.25E-02 6.68E-03 Peak BrAC 5.76E-02 1.54E-02 4.51E-02 5.66E-02 6.03E-02 2.42E-02 Peak Time 3.40E-01 7.37E-01 1.01E+00 1.52E+00 1.98E+00 5.83E-01 Area 1.19E-01 1.80E-02 1.12E-01 1.85E-01 2.85E-01 2.62E-02 Drinking Episode 7 8 9 10 11 q1 4.30E+02 5.56E+04 6.00E+02 3.40E+02 6.70E-01 q2 1.19E+00 7.70E+01 6.95E-01 1.18E+00 1.36E-02 q3 1.16E+00 5.67E+01 6.31E-01 1.16E+00 7.96E-01 q4 7.59E+03 4.03E+05 1.19E+02 8.49E+01 8.39E+06 q5 3.90E+04 2.17E+06 2.43E+02 3.76E+02 9.51E+06 TAC L 2 4.06E-02 6.98E-02 6.47E-02 3.07E-02 8.50E-02 TAC L 1 5.63E-03 7.44E-03 1.58E-02 4.53E-03 1.25E-02 BrAC L 2 2.93E-02 6.75E-02 4.11E-02 2.42E-02 7.18E-02 BrAC L 1 4.28E-03 1.63E-02 8.42E-03 4.49E-03 1.32E-02 Peak BrAC 1.62E-02 1.73E-02 2.96E-02 1.75E-02 6.12E-02 Peak Time 1.58E+00 2.67E-02 1.60E+00 1.07E+00 9.17E-01 Area 2.76E-02 1.45E-02 7.32E-02 2.49E-02 1.19E-01 Table 10.4.1: Values of the model parameters and various error values when calibrating on the 11 dierent drinking episodes. 69 A summary of the t parameters is given in Table 10.4.2. summary minimum maximum mu sigma q 1 5.59E-01 5.56E+04 5.26E+03 15934.05142 q 2 2.77E-04 7.70E+01 7.59E+00 21.95723919 q 3 6.31E-01 5.67E+01 5.95E+00 16.04691345 q 4 1.39E-01 8.39E+06 8.56E+05 2390643.526 q 5 2.01E-02 9.51E+06 1.16E+06 2718624.267 TAC L 2 2.35% 11.47% 6.16% 2.69% TAC L 1 0.45% 2.10% 1.09% 0.47% BrAC L 2 2.42% 13.12% 5.72% 3.08% BrAC L 1 0.43% 2.25% 1.19% 0.56% Peak BrAC 1.54% 6.12% 3.65% 1.88% Peak Time 0.03 1.98 1.03 0.57 Area 0.01 0.29 0.09 0.08 Table 10.4.2: Summary of data presented in Table 10.4.1 Figure 10.4.1 and Figure 10.4.2 are two plots which compare the experi- mental data to the output of the t model, that is the model with the optimal parameters, for drinking episode 1 and drinking episode 4, Figure 10.4.1: Estimated and Experimental TAC and BrAC using parameters estimated from Episode 1 70 Figure 10.4.2: Estimated and Experimental TAC and BrAC using parameters estimated from Episode 4 As the standard deviations of parameters displayed in the above table show, optimal parameters from dierent episodes are relatively far away from each other, which, at least in theory, should not be the case since these episodes are from a single subject. Therefore, we suspect that that optimal parameters are not unique This implies that the alcohol system is to some extent under- determined. To remedy this, we apply a common technique to deal with an under-determined system, a technique known as Tikhonov regularization. This will be discussed below. 10.5 Model Parameters Using Multiple Drink- ing Episodes Suppose there are m drinking episodes in the experimental data, and that training data of the i-th drinking episode isf^ y ij ; ^ u ij g ni j=1 , where ^ y ij is the j-th 71 TAC measurement, ^ u ij is the j-th measurement of amount of alcohol ingestion. Suppose further that we now use the data from several drinking episodes to estimate the model parameters q. The optimization problem we now need to solve is given by ^ q = arg min q2Q J N (q) = arg min q2Q m X i=1 ni X j=1 jy N i;j (q) ^ y ij j 2 subject to the abstract discrete system below ' N i;j+1 =S N (;q)' N i;j + ^ B N (q)F (j;' N i;j ; ^ u i;j ;q) ' N i;0 = 0 y N i;j (q) =C N (q)' N i;j ;j = 1; 2; ;n i ;i = 1; 2;:::;m; where S N (t;q) and ^ B N (q) are computed numerically in Section 10.2. This optimization problem is solved in the same manner as was the optimization problem in the previous section. The optimal parameters estimated using all M = 11 drinking episodes si- multaneously are q = [549:70; 0:1963; 0:6425; 549:87; 1010:69] Table 10.5.1 shows the dierences between the tted model produced BrAC and TAC and the experimentally observed BrAC and TAC. A summary of the above metrics is given in Table 10.5.2. We have also included a number of plots (Figure 10.5.1, Figure 10.5.2, Figure 10.5.3 and Figure 10.5.4) comparing the experimental BrAC and TAC curves to the BrAC and TAC curves produced by the t model for a number of drinking episodes. 72 Drinking Episode 1 2 3 4 5 6 7 8 9 10 11 TAC L 2 0.0472 0.0575 0.0767 0.1700 0.1099 0.0531 0.0500 0.0358 0.1115 0.0627 0.1489 TAC L 1 0.0110 0.0062 0.0126 0.0209 0.0198 0.0065 0.0053 0.0056 0.0167 0.0093 0.0221 BrAC L 2 0.0916 0.0788 0.1185 0.1760 0.1479 0.1045 0.0935 0.0821 0.1544 0.0615 0.1339 BrAC L 1 0.0199 0.0120 0.0198 0.0198 0.0202 0.0147 0.0096 0.0128 0.0190 0.0065 0.0165 Peak BrAC 0.0445 0.0154 0.0389 0.0396 0.0569 0.0088 0.0097 0.0096 0.0395 0.0092 0.0485 Peak Time 0.34 0.78 1.4567 1.5167 1.98 0.5833 1.58 0.8533 1.9733 1.0933 0.9167 Area 0.1394 0.0527 0.1560 0.1841 0.3095 0.0258 0.0304 0.0303 0.1635 0.0302 0.1756 Table 10.5.1: Various error values when calibrating on the 11 dierent drinking episodes summary minimum maximum mu sigma TAC L 2 0.035774667 0.170086756 0.083964973 0.042673141 TAC L 1 0.005309634 0.022119586 0.012365812 0.006194656 BrAC L 2 0.061469461 0.176004853 0.112983686 0.034437465 BrAC L 1 0.00654322 0.020180802 0.015533849 0.00452998 Peak BrAC 0.008778809 0.056932697 0.029144072 0.017729151 Peak Time 0.34 1.98 1.188484848 0.524287614 Area 0.025760715 0.309539765 0.117964174 0.087477336 Table 10.5.2: Summary of data presented in Table 10.5.1 Figure 10.5.1: Estimated and experimental TAC and BrAC for drinking episode 1 using parameters estimated using drinking episodes 1-7 73 Figure 10.5.2: Estimated and experimental TAC and BrAC for drinking episode 4 using parameters estimated using drinking episodes 1-7 Figure 10.5.3: Estimated and experimental TAC and BrAC for drinking episode 5 using parameters estimated using drinking episodes 1-7 74 Figure 10.5.4: Estimated and experimental TAC and BrAC for drinking episode 11 using parameters estimated using drinking episodes 1-7 10.6 Estimating Model Parameters with Regu- larization Even though the primary focus of this thesis is to investigate the quantitative characterization of drinking patterns and behavior (i.e. the drinking parame- ters), rather than estimating the model parameters, it would be desirable if the model parameters used when the drinking parameters are estimated exhibit some degree or level of uniqueness and identiability. This can be achieved by introducing Tikhonov regularization into the objective function when the model parameters are estimated. We use the same training data as in the previous sec- tion, but now we add a weighted regularization term to the objective function 75 as in, ^ q = arg min q2Q J N (q) = arg min q2Q m X i=1 ni X j=1 jy N ij (q) ^ y ij j 2 +wjjqq 0 jj 2 where ! is the regularization weight or coecient and q 0 is an initial guess, or nominal values for the parameters based on any available prior knowledge. Numerically this optimization problem is solved in the same manner as those in the previous section. The adjoint method for computing the gradient is easily modied to account for the added regularization term. Below is a summary of the parameters estimated from each drinking episode along with plots comparing experimental to modeled (using the parameters t with regularization present) BrAC and TAC curves for a number of episodes. Drinking Episode q 1 q 2 K m 1 0.91 0.64 0.81 0.62 1.39 2 0.98 1.21 0.97 1.22 1.28 3 1.05 0.20 0.81 0.24 0.22 4 0.59 0.46 0.83 0.56 1.00 5 1.06 0.11 0.74 0.14 0.02 6 0.97 0.86 1.28 0.88 1.17 7 0.98 1.07 1.34 1.06 1.12 8 0.97 1.17 1.46 1.15 1.12 9 0.75 0.67 0.60 0.69 1.35 10 0.97 1.16 1.36 1.19 1.33 11 0.79 0.67 0.75 0.63 1.40 Table 10.6.1: Values of the model parameters when calibrating on the 11 dier- ent drinking episodes 76 Figure 10.6.1: Estimated and experimental TAC and BrAC using parameters estimated based on drinking episode 1 data with regularization. Figure 10.6.2: Estimated and experimental TAC and BrAC using parameters estimated based on drinking episode 7 data with regularization. Using the rst 7 drinking episodes as data, and now including regularization, 77 the optimal parameters are given by q 1 = 0:99;q 2 = 0:58;b = 0:76;K = 0:90;m = 2:09: Note that we will use this set of parameters in our drinking pattern and behavior (i.e. estimation of drinking parameters) study in a later chapter of the thesis. Below are several plots to comparing experimental BrAC and TAC curves to BrAC and TAC curves produced by the t model with the parameters t using regularization. Figure 10.6.3: Estimated and experimental TAC and BrAC for drinking episode 1 using parameters estimated using drinking episodes 1-7 but now with regular- ization 78 Figure 10.6.4: Estimated and experimental TAC and BrAC for drinking episode 5 using parameters estimated using drinking episodes 1-7 but now with regular- ization Figure 10.6.5: Estimated and experimental TAC and BrAC for drinking episode 11 using parameters estimated using drinking episodes 1-7 but now with regu- larization 79 Chapter 11 Estimation of Random Dynamical Systems Our general approach is based on two relatively recent treatments by Banks and his co-authors (see [2, 5]) in which the authors developed an abstract ap- proximation framework for the estimation of random dynamical systems. We begin by summarizing the results from those treatments that we will require in the sequel and then we will apply those results to the particular problem and model of interest to us here. 11.1 The Prohorov Metric on a Set of Probabil- ity Measures Letf;dg be a compact Hausdor metric space, let C B () be the space of bounded, continuous, and not necessarily linear real valued functionals on , and let C B () denote its (continuous) dual space. Let P() denote the set of probability measures on a sigma algebra of measurable subsets off;dg. A theorem of Riesz guarantees that for all f 2 C B () there exists a unique nite signed Borel measure onf;dg such that f (f) = Z f()d(); for all f2C B (), withjjf jj =jj(). It follows from this result that we may identifyP() with a subset of C B () , via the mapping :P()!C B () 80 where for P2P(), (P ) =f P 2C B () , withjj (P )jj =jjf P jj =P () = 1, (P )(f) =f P (f) = Z f()dP (); for f2 C B (), and where the range of is given byR( ) =ff 2 C B () : f (f) 0;f2C B ();jjf jj = 1g. Denition 11.1.1. ForfP n g;P2P() we say the sequencefP n g converges weakly to P and write P n wk* !P if f Pn converges pointwise to f P . We note that even though the denition above is for weak convergence of probability measures, the use of the notation P n wk* ! P is justied because point wise convergence in C B () will follow if f Pn converges to f P in the weak* topology. Indeed, if we let denote the canonical map from C B () into C B () , and for eachf 2C B () ,f (f Pn )!f (f P ), then it follows that for each f2C B () f Pn (f) =(f)(f Pn )!(f)(f P ) =f P (f): We note that if C B () is re exive, then is surjective and f Pn wk* !f P would be equivalent to f Pn converging pointwise to f P . The local neighborhood system dened by B " (P ) =fQ2P() :jf Q (f)f P (f)j<";f2C B ()g where P2P(), and " > 0 can be shown to generate the weak* topology on P(). Moreover, this topology is metrizable. Denition 11.1.2. ForP;Q2P(), the Prohorov metric,, onP() is given 81 by (P;Q) = inff"> 0 :Q(A)P (A " )+" and P (A)Q(A " )+"; for all A2 g; where P;Q2P() and for all " > 0, and A2 , A " is an "-neighborhood of A dened by A " =f2 : d(A;) < "g, if A is nonempty, and A " =; if A =;. The following theorem is proven in [2] and [5]. Theorem 11.1.1. Iff;dg is a separable metric space, then is a metric on P(). Moreover, iffP n gP() and P2P(), then P n wk* !P as n!1 if and only if (P n ;P )!P as n!1, and consequently, the topology induced by the metric is the same as the weak topology of measures. The following key properties and consequences of the Prohorov metric play an important role in applying the Prohorov framework to the estimation problem of interest to us here. The proofs of these results can be found in [2, 5]. Denition 11.1.3. For2 the Dirac measure or distribution, , on is given by (A) = A (), forA2 , where A denotes the standard character- istic function of the set A. LetDP() denote the set of all Dirac measures, D =f :2 g. Proposition 11.1.1. Iff;dg is a separable metric space( ; ) = minfd(;); 1g. Proposition 11.1.2. Iff;dg is a separable metric spacef n g Cauchy in f;dg if and only iff n g Cauchy infP();g. Proposition 11.1.3. Iff;dg is a separable metric space the setD is sequen- tially closed infP();g. Proposition 11.1.4. Iff;dg is a separable metric space, it is complete if and only iffP();g is a complete metric space. 82 Proposition 11.1.5. Iff;dg is a separable metric space, it is compact if and only iffP();g is a compact metric space. Proposition 11.1.6. Iff;dg is a separable metric space, the metric space fP();g is also separable, and moreover, if S =f k g 1 k=1 is dense inf;dg, then the set fP2P() :P = K X k=1 p k k ; k 2S; k 2D;K2N;p k 2 [0; 1]\Q; K X k=1 p k = 1g is countable and dense infP();g. 11.2 Estimation of Random Discrete Time Dy- namical Systems In this section, we adapt the results in [2] and [5] on the estimation of random dynamical systems using the Prohorov metric outlined in the previous section to discrete or sampled time innite dimensional systems. We combine this with some recent results in [29, 30, 31, 32] on the estimation of random parameters in abstract parabolic systems. We consider the family of discrete or sampled time initial value problems that are set in an, in general, innite dimensional Hilbert state space,H, given by x j+1;i =g(t j ;x j;i ;u i ;u); j = 0;:::;n i ; i = 1; 2;:::;m; (11.2.1) x 0;i =x 0;i ; i = 1; 2;:::;m; (11.2.2) whereg :R + H Q ni j=0 R R !H and forj = 0;:::;n i andi = 1; 2;:::;m, u i =fu i;j g is an external input or control with u i;j 2 R , and t j = j, with > 0 the length of the sampling interval, describing the dynamics of a process 83 common to the entire population. In equations (11.2.1) and (11.2.2) the index i refers to an individual subject, while j denotes time. The parameter vector u which takes values in a compact subset, U of R determines how the input to the system, u i;j 2 R , and t j = j, aects or drives the dynamics. These are the quantities or parameters (or simply, a parameter) that we want to estimate based on observed data. In our treatment below we will assume that u is random and rather than estimate it directly, it is its distribution in the form of a joint probability density function, or joint cumulative distribution function, or probability measure, that will be estimated. We assume that we can observe some function of the solutions of (11.2.1)- (11.2.2), x j;i , as given by the output equation y j;i =y(t j ;x 0;i ;u i ;u) =h(x j;i ;x 0;i ;u i ;u); j = 0;:::;n i ; i = 1; 2;:::;m; (11.2.3) where h :HH Q ni j=0 R R !R . In equations (11.2.1)-(11.2.3), we assume u is in some compact subset U of R , the set of admissible values foru. We assume that the values of these param- eters are specic to each individual in the population. Therefore, assuming that the parameters,u, are samples from a random vectorU, dened on a probability spacef; ;P 0 g, the objective is to estimate their (joint) distribution based on the data sampled from the population. For this purpose, we assume that the distribution of these random vectors is described by a joint cdf F 0 2F, where F represents a set of feasible cdfs. We consider two approaches to estimating F 0 , 1) by estimating the probability measure, P 0 , in the underlying probability space and 2) by assuming F 0 is absolutely continuous and estimating the joint probability density function. We consider two possible statistical models, the naive pooled data model and 84 the mixed eects model (see, for example, [15, 16, 17, 31]) wherein randomness in the parameters, u, are used to quantify uncertainty between subjects, and randomness in the output or measurements, y j;i given in (11.2.3) is intended to capture uncertainty within individual subjects. In the case of the mixed eects model, we assume that the observed data points are of the form V j;i =y j;i +" j;i ; j = 0;:::;n i ; i = 1;:::;m; where " j;i ; j = 0;:::;n i ; i = 1;:::;m, representing measurement noise are as- sumed to be independent across subjects (i.e with respect to i), conditionally independent with respect to U within subjects (i.e with respect to j), iden- tically distributed with mean 0 and known common variance 2 , and with " j;i '; j = 0;:::;n i ; i = 1;:::;m. Letting V =fV j;i g, then V can be thought of as a vector inq-dimensional Euclidean space, R q , for someq, and under these assumptions, the likelihood takes the form L(F ;V ) = m Y i=1 Z 1 1 L i (u;V )dF (u) = m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x j;i ;x 0;i ;u i ;u))dF (u); (11.2.4) or L(P ;V ) = m Y i=1 Z L i (U();V )dP () = m Y i=1 Z ni Y j=0 '(V j;i h(x j;i ;x 0;i ;u i ;U()))dP (); (11.2.5) whereF (u) =P (fUug). One could then seek a maximum likelihood estima- tor forF 0 by maximizing either (11.2.4) overF2F, or (11.2.5) overP2P(). Solving the optimization problems for the MLE presents a number of signif- 85 icant computational challenges including the ecient evaluation of an integral over a high dimensional set, under ow issues when dealing with products of a large quantity of small numbers, etc. Of principal concern to us here however, is the discretization of the innite dimensional state and feasible set and as- sociated convergence questions. When the MLE problem is solved, one would actually be maximizing the functionals L N;M (F M ;V ) = m Y i=1 Z 1 1 L N i (u;V )dF M (u) = m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;u))dF M (u); (11.2.6) or L N;M (P M ;V ) = m Y i=1 Z L N i (U();V )dP M () = m Y i=1 Z ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;U()))dP M (); (11.2.7) where the positive integers N and M denote discretization indices. Once one uses (11.2.6) or (11.2.7) to respectively nd approximating MLE's, F N;M MLE or P N;M MLE , the natural question to ask is if in some senseF N;M MLE !F MLE orP N;M MLE ! P MLE , as N;M!1. One could also ask questions about the consistency of these estimators. As is typically the case in MLE problems, one would make use of the log function so as to transform some of the products in (11.2.6) or (11.2.7) into sums and then under appropriate regularity assumptions, a gradient based optimization scheme could be used to solve the maximization problems. An E- M approach might also work wherein observations ofU are treated as missing data (see, for example, [11]). Alternatively, we could assume that we have aggregate data and seek the naive pooled data estimator. In this case, we assume thatn i =n,x 0;i =x 0 and 86 u i = u;i = 1; 2; ;m, a statistical model of the form V j =E[y j (U)jP 0 ] +" j ; j = 0;:::;n (11.2.8) where in (11.2.8)," j ; j = 0;:::;n, represent measurement noise and are assumed to be independent and identically distributed with mean 0 and known common variance 2 . Once again lettingV =fV j g, thenV can be thought of as a vector inq-dimensional Euclidean space, R q , for someq, and forF2F, orP2P() dene v(t j ;F ) =E[y(t j ;x 0 ; u;U)jjF ] = Z 1 1 h(x j ;x 0 ; u;u)dF (u); (11.2.9) or v(t j ;P ) =E[y(t j ;x 0 ; u;U)jjP ] = Z h(x j ;x 0 ; u;U())dP (); (11.2.10) the mean behavior at time t j , j = 0;:::;n, ifUF . The estimation problem is to then estimate the cdf,F 0 , using a least squares approach with either ^ F = arg min F2F J(F ;V ) = arg min F2F n X j=0 (V j v(t j ;F )) 2 : (11.2.11) or ^ P = arg min P2P() J(P ;V ) = arg min P2P() n X j=0 (V j v(t j ;P )) 2 : (11.2.12) where the v(t j ;F ) are as given in (11.2.9) and the v(t j ;P ) are as given in (11.2.10). Once again, solving the optimization problems given in either (11.2.11) or 87 (11.2.12) will typically require nite dimensional approximation of the dynam- ical system given in (11.2.1)-(11.2.2), and the parameterization of the feasible set of cdfs,F or measuresP(). As in the MLE case, the problems we actually end up solving are given by v N;M (t j ;F M ) =E[y(t j ;x N 0 ; u;U)jjF M ] = Z 1 1 h(x N j ;x N 0 ; u;u)dF M (u); (11.2.13) or v N;M (t j ;P M ) =E[y(t j ;x N 0 ; u;U)jjP M ] = Z h(x N j ;x N 0 ; u;U())dP M (); (11.2.14) the mean behavior at time t j , j = 0;:::;n, ifUF M . The estimation problem is to then estimate the cdf,F 0 , using a least squares approach with either ^ F N;M = arg min F M 2F J N;M (F M ;V ) = arg min F M 2F n X j=0 (V j v N;M (t j ;F M )) 2 : (11.2.15) or ^ P N;M = arg min P M 2P() J N;M (P M ;V ) = arg min P M 2P() ni X j=0 (V j v N;M (t j ;P M )) 2 : (11.2.16) where once again the positive integers N and M denote discretization indices, v N;M (t j ;F M ) are as given in (11.2.13) and the v N;M (t j ;P M ) are as given in (11.2.14). As in the MLE case, the questions that one might ask are: under what assumptions will it follow that ^ F N;M ! ^ F , ^ P N;M ! ^ P , as N;M!1, and are these estimators consistent. 88 If we are willing to restrict our search to absolutely continuous distributions, both the mixed eects and naive pooled data estimation problems can be re- formulated in terms of densities. In the mixed eect case, (11.2.4) and (11.2.6) become L(f;V ) = m Y i=1 Z 1 1 L i (u;V )f(u)du = m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x j;i ;x 0;i ;u i ;u))f(u)du; (11.2.17) and L N;M (f M ;V ) = m Y i=1 Z 1 1 L N i (u;V )f M (u)du = m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;u))f M (u)du; (11.2.18) respectively, and in the naive pooled data case, (11.2.9), (11.2.11), (11.2.13), and (11.2.15) become v(t j ;f) =E[y(t j ;x 0 ; u;U)jjf] = Z 1 1 h(x j ;x 0 ; u;u)f(u)du; (11.2.19) ^ f = arg min f2f J(f;V ) = arg min f2f n X j=0 (V j v(t j ;f)) 2 ; (11.2.20) v N;M (t j ;f M ) =E[y(t j ;x N 0 ; u;U)jjf M ] = Z 1 1 h(x N j ;x N 0 ; u;u)f M (u)du; (11.2.21) 89 and ^ f N;M = arg min f M 2f J N;M (f M ;V ) = arg min f M 2f n X j=0 (V j v N;M (t j ;f M )) 2 ; (11.2.22) respectively. In these cases the convergence questions become f N;M MLE ! f MLE and ^ f N;M ! ^ f, as N;M !1. We note that the pdf estimation problems can be approached either parametrically (e.g. estimating the mean, , and standard deviation, , of a normal distribution, or the mean, , of a Poisson distribution), or non-parametrically, where it is assumed, for example that f M is a linear combination of spline functions (with appropriate constraints) and it is the coecients in the linear combination that are estimated. In the lat- ter case, it is eectively the shape of f 0 that is being estimated. In the next two sections we summarize the existence, consistency and convergence results for non-parametrically estimating the underlying probability measure and para- metrically estimating the probability density function. 11.3 Using the Prohorov Metric Framework to Non-ParametricallyEstimatetheUnderly- ing Probability Measure In this section we present three theorems without proof from [2] and [5] that summarize their existence, consistency, and convergence results for the case where we are estimating the probability measure. The proofs can be found in those papers. Theorem 11.3.1. (Existence of the Estimators) Letf;dg be a separable and compact metric space and assume that for each P 2fP();g the mappings given by (11.2.4) and (11.2.13), L(P ;)R q ! R and J(P ;) : R q ! R are 90 measurable and that for each V 2 R q the mappings again given by (11.2.4) and (11.2.13), L(;V ) : fP();g ! R and J(;V ) : fP();g ! R are continuous. Then there exists a measurable functions P MLE : R q !P() and ^ P : R q !P() such thatL(P MLE (V );V ) = sup P2P() L(P ;V ), for V 2 R q , and J( ^ P (V );V ) = inf P2P() J(P ;V ), for V 2R q . We note that an analogous result holds for the nite dimensional estima- tion problems dened in terms of (11.2.7) and (11.2.16). Not surprisingly, the essence of the proof is maximization or minimization of continuous functions over compact sets and the various properties of the Prohorov metric onP(). In [2] and [5] the authors were also able to establish a consistency result for the naive pooled data estimator. They require the following additional assump- tions: (a) The measurement noisef" j g is i.i.d. with respect to a probability space f ; ;P g with E[" j jjP ] = 0 and Var[" j jjP ] = 2 , (b) The metric spacef;dg is separable and compact, and we consider the space of probability measures on ,P() to be endowed with the Prohorov metric . (c) Fori = 1; 2;:::,n i =n andn =T for some positive integer n and some T > 0, where is the sampling time dened earlier (d) For each P2P(), The modelfv(t j ;P )g given by (11.2.10) converges to continuous functions v(;P ) on [0;T ], i = 1; 2;:::;m as the sampling interval ! 0 with v continuous fromfP();g into C[0;T ]. (e) There exists a measure on [0;T ] such that for all g2C[0;T ] 1 n n X j=1 g(t j ) = Z T 0 g(t)d n (t)! Z T 0 g(t)d(t);n!1: 91 (f) P 0 2P() is the unique minimizer of J 0 inP() where J 0 (P ) = 2 + Z T 0 (v(t;P 0 )v(t;P )) 2 d(t): (11.3.1) Then a straight forward application of Theorem 4.2 in [5] can then be used to establish the following lemma and theorem (see [31]). Lemma 11.3.2. If Assumptions (a) - (f) above hold, then there exists an event A2 with P (A) = 1 such that for all !2A and J as given in (11.2.12) we have 1 n J(P ;V )!J 0 (P ); as n!1 and ! 0, with n =T , uniformly in P for P2P(), where J is given by J(P ) = n X j=0 (V j v(t j ;P )) 2 ; and J 0 by (11.3.1). Theorem 11.3.3. (Consistency of the Estimator ^ P ) Let ^ P2P() be as de- ned in (11.2.12). Then under the assumptions of the previous Lemma the estimator ^ P2P() is consistent for P 0 . That is ^ P! P 0 in probability with respect to the probability measure P , as m;n!1, and ! 0 with n =T . Theorem 11.3.4. (Convergence of the Approximating Estimators) Letf;dg be a separable and compact metric space and let V 2 R q be given and xed. Assume (1) that the mappings given by (11.2.7) and (11.2.16), L N;M (;V ) : fP();g! R and J N;M (;V ) :fP();g! R are continuous; (2) For any sequencefP M g2fP();g withP M !P infP();g, we havew N;M i (P M )! w i (P ) as N;M !1, i = 1; 2;:::;m, where, recalling (11.2.7) and (11.2.5), 92 w N;M i (P M ) and w i (P ) are given by w N;M i (P M ) = Z L N i (U();V )dP M (); and w i (P ) = Z L i (U();V )dP (); fori = 1; 2;:::;m, respectively, and we havev N;M (t j ;P M )!v(t j ;P ) asN;M! 1,j = 0; 1; 2;:::;n, wherev N;M (t j ;P M ) andv(t j ;P ) are given by (11.2.10) and (11.2.14), respectively; and (3)w i (P ),i = 1; 2;:::;m,v(t j ;P ),j = 1; 2;:::;n, are uniformly bounded for all P2fP();g. Then, there exists optimizers, P N;M MLE and ^ P N;M for the optimization problems given in (11.2.7) and (11.2.16), re- spectively. In addition, the sequencesfP N;M MLE g andf ^ P N;M g admit subsequences which converge to optimizers for the optimization problems given in (11.2.5) and (11.2.12), respectively. It is noted in [2] and [5] that their Prohorov metric framework for estimating the distribution of random parameters in dynamical systems is computationally constructive. The space of probability measures P() can be discretized by considering measures that are linear combinations of Dirac measures of the form P M (A) = P M k=1 p M k M k (A), for A2 where the nodesf M k g M k=1 and the weightsfp M k g M k=1 R + with P M k=1 p M k = 1. In the case of the mixed eects model, the likelihood to be maximized takes the form 93 L N;M (P M ;V ) = m Y i=1 Z L N i (U();V )dP M () = m Y i=1 Z ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;U()))dP M () = m Y i=1 Z ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;U()))d M X k=1 p M k M k () = m Y i=1 M X k=1 p M k ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;U( M k ))): (11.3.2) In the case of the naive pooled data estimator, we have that v N;M (t j ;P M ) =E[y(t j ;x N 0 ; u;U)jjP M ] = Z h(x N j ;x N 0 ; u;U())dP M () = Z h(x N j ;x N 0 ; u;U())d M X k=1 p M k M k () = M X k=1 p M k h(x N j ;x N 0 ; u;U( M k )); (11.3.3) in which case J N;M (P M ;V ) = n X j=0 (V j v N;M (t j ;P M )) 2 = n X j=0 (V j M X k=1 p M k h(x N j ;x N 0 ; u;U( M k ))) 2 : (11.3.4) The optimization can be with respect just thefp M k g M k=1 with thef M k g M k=1 chosen a-priori and xed, or the optimization can be with respect to both the fp M k g M k=1 and thef M k g M k=1 . In the case where the statistical model is the naive pooled data estimator and the optimization is over only thefp M k g M k=1 , the op- 94 timization problem given in (11.3.3) and (11.3.4) becomes one of minimizing a quadratic form over a subset of M dimensional Euclidean space, ^ P N;M = arg min P M 2P() J N;M (P M ;V ) = arg min p M 2R + M; P M k=1 p M k =1 f(p M ) T Hp M +L T p M g (11.3.5) where H k;l = n X j=0 h(x N j ;x N 0 ; u;U( M k ))h(x N j ;x N 0 ; u;U( M l ));k;l = 1; 2;:::;M; (11.3.6) and L k =2 n X j=0 V j h(x N j ;x N 0 ; u;U( M k ));k = 1; 2;:::;M: (11.3.7) With regard to verifying the hypotheses of the convergence theorem, The- orem (11.3.4), so long as g : R + H Q ni j=0 R R ! H, h : HH Q ni j=0 R R ! R , andU : ! R are continuous and since we assumed that is compact and therefore thatfP;g is compact, conditions (1) and (3) are easily veried. It is shown in [2] and [5] that in order to verify condition (2), it need only be argued that y N j;i =y(t j ;x N 0;i ;u i ;U()) =h(x N j;i ;x N 0;i ;u i ;U())!y j;i ;N!1; uniformly in for 2 , and that the mappings ! g(t j ;x j;i ;u i ;U()) and !h(x N j;i ;x N 0;i ;u i ;U()) are continuous on . We note that in actual practice, it will often be more convenient to choose the sample space, as the range ofU,R(U), rather than as the domain ofU,D(U), and to then estimate a measure P on = R(U) . In this caseU() =, for 2 and we would let P M (A) = P M k=1 p M k u M k (A), for A2 R(U) where the 95 nodesfu M k g M k=1 R(U) and the weightsfp M k g M k=1 R + with P M k=1 p M k = 1. Then (11.3.2) - (11.3.7) become (11.3.8) - (11.3.13) L N;M (P M ;V ) = m Y i=1 M X k=1 p M k ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;u M k )); (11.3.8) v N;M (t j ;P M ) = M X k=1 p M k h(x N j ;x N 0 ; u;u M k ); (11.3.9) in which case J N;M (P M ;V ) = n X j=0 (V j M X k=1 p M k h(x N j ;x N 0 ; u;u M k )) 2 ; (11.3.10) and as before ^ P N;M = arg min P M 2P(R(U)) J N;M (P M ;V ) = arg min p M 2R + M; P M k=1 p M k =1 f(p M ) T Hp M +L T p M g; (11.3.11) but now with H k;l = n X j=0 h(x N j ;x N 0 ; u;u M k )h(x N j ;x N 0 ; u;u M l );k;l = 1; 2;:::;M; (11.3.12) and L k =2 n X j=0 V j h(x N j ;x N 0 ; u;u M k );k = 1; 2;:::;M: (11.3.13) 96 11.4 Parametrically Estimating the Probability Density Function In this section we discuss results for parametrically estimating the distribu- tion in the form of a probability density function. These results are essentially the same as those developed in [29, 30, 31, 32]. Let denote the set of fea- sible parameters, a suset of Euclidean r-space, R r and typically assumed to be compact. Let f denote the set of feasible probability density functions, f parameterized by 2 . Recalling (11.2.17), (11.2.18), (11.2.19), (11.2.20), (11.2.21), and (11.2.22) we then consider the following optimization problems over and their approximations: ^ = arg max 2 L(f ;V ) = arg max 2 m Y i=1 Z 1 1 L i (u;V )f (u)du = arg max 2 m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x j;i ;x 0;i ;u i ;u))f (u)du; (11.4.1) ^ N = arg max 2 L N (f ;V ) = arg max 2 m Y i=1 Z 1 1 L N i (u;V )f (u)du = arg max 2 m Y i=1 Z 1 1 ni Y j=0 '(V j;i h(x N j;i ;x N 0;i ;u i ;u))f (u)du; (11.4.2) 97 ^ = arg min 2 J(f ;V ) = arg min 2 n X j=0 (V j v(t j ;f )) 2 ; (11.4.3) ^ N = arg min 2 J N (f ;V ) = arg min 2 n X j=0 (V j v N (t j ;f )) 2 : (11.4.4) We then have the following existence and convergence result for the above de- ned estimators. Theorem 11.4.1. Let R r be compact. If A. The maps on ,7!f (u), for almost everyu2R , and7!J N (f ;V ), for all N and f 2f are continuous, B. For any sequence of densities f M 2f with lim M!1 f M (u) = f (u), a.e. u2 R , for some f 2f , we have v N (t j ;f M ) converging to v(t j ;f ) for all j2f0;:::;ng as N;M!1, and C. The v(t j ;f ) and v N (t j ;f ) are uniformly bounded for all j2f0;:::;ng and f 2f , then it will follow that there exist solutions ^ N to the estimation problems over , given in either (11.4.2) or (11.4.4), and there exists a subsequence of the ^ N 's that converges to a solution ^ of the estimation problems over given in either (11.4.1) or (11.4.3), respectively. As in the previous section, With regard to verifying the hypotheses of the convergence theorem, Theorem (11.4.1), so long asg :R + H Q ni j=0 R R ! H andh :HH Q ni j=0 R R !R are continuous and since we assumed that andU are compact and therefore thatf is compact, conditions (A) and (C) are easily veried. In order to verify condition (B), it need only be argued 98 that y N j;i =y(t j ;x N 0;i ;u i ;u) =h(x N j;i ;x N 0;i ;u i ;u)!y j;i ;N!1; for each u2 U, i = 1; 2;:::;m, j = 0; 1; 2;:::;n i , where y j;i is given in (11.2.3), and that the mappings u!g(t j ;x j;i ;u i ;u) and u!h(x N j;i ;x N 0;i ;u i ;u) are con- tinuous on U. It is also possible to establish a consistency result for the estimator ^ 2 . We require the following additional assumptions: (a) The measurement noisef" j g is i.i.d. with respect to a probability space f ; ;P g with E[" j jjP ] = 0 and Var[" j jjP ] = 2 , (b) The feasible set of parameters is compact (i.e. closed and bounded since it is nite dimensional) and has nonempty interior, (c) Fori = 1; 2;:::;n andn =T for some positive integern and someT > 0, where is the sampling time, (d) ThatV j =E[y j jjf 0 ]+" j , for some 0 2intfg, where forj = 0; 1; 2;:::;n y j is given by (11.2.3) (e) For each f 2f , The modelfv(t j ;f )g given by (11.2.19) converges to continuous functions v(;f ) on [0;T ] as the sampling interval ! 0 with the map !v(;f ) continuous from into C[0;T ]. (f) There exists a measure on [0;T ] such that for all g2C[0;T ] 1 n n X j=1 g(t j ) = Z T 0 g(t)d n (t)! Z T 0 g(t)d(t);n!1: (f) The density f 0 2 is the unique minimizer of J 0 inf where J 0 (f ) = 2 + Z T 0 (v(t;f 0 )v(t;f )) 2 d(t); (11.4.5) Then a straight forward application of Theorem 4.2 in [5] can then be used to establish the following lemma and theorem (see [31]). 99 Lemma 11.4.2. If in addition to Assumptions (a) { (e) above we assume that the maps 7! f (u) from to R are continuous for a.e. u2 U, then there exists an event A2 with P (A) = 1 such that for all !2A and J as given in (11.2.20) we have 1 n J(f ;V )!J 0 (f ); as n!1 and ! 0, with n = T , uniformly in for 2 , where J 0 by (11.4.5) and J is given by (11.4.6) J(f ;V ) = n X j=0 (V j v(t j ;f )) 2 ; (11.4.6) Theorem 11.4.3. (Consistency of the estimator ^ Let ^ 2 be as dened in (11.4.3). Then under the assumptions of Lemma (11.4.2), the estimator ^ 2 is consistent for 0 . That is ^ rightarrow 0 in probability with respect to the probability measure P , as n!1, and ! 0 with n =T . 100 Chapter 12 Estimation of Random Parameters in a Distributed Parameter Model for the Metabolism and Transport of Ethanol in the Human Body Based on Transdermal Alcohol Concentration In this chapter we will use the semilinear alcohol transport and metabolism model developed in the previous chapters to quantitatively characterize drinking patterns and behaviours. In order to adapt the theoretical results we introduced in the last chapter, we want to construct a random discrete time dynamical system and formulate the drinking behavior estimation problem as a problem of identifying random parameters in such system. We begin by looking at the discrete time alcohol metabolism and transport model that we developed in Chapter 6 and consider the following discrete-time initial value problem with sampling time given by ' j+1;i =S(;q)' j;i + ^ B(q)F (i;' j;i ;u j;i ;u i ;q) ' 0;i = 0 y j;i =C(q)' j;i ;j = 1; 2; ;n i ;i = 1; 2; ;m; (12.0.1) where S(;q) is a semigroup generated by A(q) dened in (4.3.2) and ^ B(q) = R 0 S()d. In the system (12.0.1), the index i refers to an individual subject, which in our case here, refers to a drinking episode, while j denotes time index (t j = 101 j). We assume that the model parameters q have already been estimated or identied using the deterministic nonlinear least squares techniques discussed earlier in the thesis and are specic to the individual subject whose drinking pattern or behavior is to be estimated. The input u j;i denotes the amount of alcohol ingested by the subject in time interval j during drinking episode i in units of standard drinks. We assume that the amount of alcohol ingested during any sampling time interval is constant over the course of the i-th drinking episode. The parameter vector u i which takes values in a compact subset of Euclidean space,U ofR p , for somep, determines the level of this constant value in the i-th drinking episode as a function of the number of drinks consumed and/or the duration of the drinking episode. So, for example if thei-th drinking episode lasts T i hours, d i standard drinks are ingested during the i-th drinking episode, and the length of the sampling interval is hours, then u i = d i =T i standard drinks per hour. If we let i = (T i ;d i ) and u j;i = , for all i;j, then the amount of alcohol ingested in each sampling interval in units of standard drinks is given by u i u i;j =u( i )u i;j = (d i =T i )u i;j = (d i =T i ) standard drinks. In what is to follow we will assume that either or bothT i andd i are random and sampled from random variablesT andd, respectively, in which case both i andu i are random and are sampled from the random vector and the random variableU =u(), respectively. In this case we have u i U() and our model now takes the form of the random dynamical system given by ' j+1;i =S(;q)' j;i + ^ B(q)F (i;' j;i ;u j;i ;U();q) ' 0;i = 0 y j;i =C(q)' j;i ;j = 1; 2; ;n i ;i = 1; 2; ;m; (12.0.2) which is now of the general form given in (11.2.1)-(11.2.3) in the previous chap- ter. 102 Now to illustrate how the estimation problem will actually be formulated in the spirit of the approach developed in the previous chapter, we might for example assume that the number of drinks ingested by the subject during a drinking episode is Poisson distributed with parameter, sod Poisson() and the duration of a drinking episode has, for example, a normal distribution with mean and standard deviation ,TN(;) or a gamma distribution,T Gamma(;), or possibly a -squared distribution with degrees of freedom, T(). We might also assume, for example, thatd andT are independent or possibly that they are in fact dependent. If we want to parametrically estimate the joint density of d and T, we might assume thatd Poisson() andT Gamma(;) and that they are independent. In this case we have = (;;)2 , denoting the set of feasible parameters, a subset of R 3 which is assumed to be compact. In this case we have f (d;T ) =f (;;) (d;T ) = d e d! 1 () T 1 e T : (12.0.3) Then in order to estimate the parameters , as in the previous chapter, we consider two possible statistical models, the naive pooled data model and the mixed eects model. In the case of the mixed eects model, we assume that the observed data points are of the form V j;i =y j;i +" j;i ; j = 0;:::;n i ; i = 1;:::;m; where" j;i ; j = 0;:::;n i ; i = 1;:::;m, representing measurement noise, identically distributed with mean 0 and known common variance 2 , and with" j;i '; j = 0;:::;n i ; i = 1;:::;m. Letting V =fV j;i g, then V can be thought of as a vector 103 inq-dimensional Euclidean space, R q , for someq, and under these assumptions, we want to solve the following optimization problem and its nite dimensional approximation, ^ = arg max 2 L(f ;V ) = arg max 2 m Y i=1 1 X d=0 Z 1 0 ni Y j=0 '(V j;i y j;i )f (d;T )dT; ^ N = arg max 2 L N (f ;V ) = arg max 2 m Y i=1 1 X d=0 Z 1 0 ni Y j=0 '(V j;i y N j;i )f (d;T )dT: On the other hand, for the naive pooled data estimator, we assume a statis- tical model of the form V j =E[y j (U)jf 0 ] +" j ; j = 0;:::;n; for some 0 2 , where " j ; j = 0;:::;n, represent measurement noise and are assumed to be independent and identically distributed with mean 0 and known common variance 2 . The data, V =fV j g, can be thought of as a vector in q-dimensional Euclidean space, R q , for some q, and for f 0 with 0 2 . Then let v(t j ;) =E[y j (U)jf ] = 1 X d=0 Z 1 0 y j f (d;T )dT; be the mean behavior of the model at time t j , j = 0;:::;n. The estimation problem is to then estimate the parameters, , using a least 104 squares approach with ^ = arg min 2 J(f ;V ) = arg min 2 n X j=0 (V j v(t j ;f )) 2 : The nite dimensional optimization problem in this case is ^ N = arg min 2 J N (f ;V ) = arg min 2 n X j=0 (V j v N (t j ;f )) 2 : The existence and consistency of the above estimators follow from Theorem 11.4.1. In the case of the Prohorov metric framework for estimating the distribution of random parameters in dynamical systems, Let =f0; 1; 2; 3;:::; ^ dg [0; ^ T ), where ^ d and ^ T are suciently large and xed (recall, we want to be compact), and then as we noted in the previous chapter, the space of probability measures P() can be discretized by considering measures that are linear combinations of Dirac measures of the form P M (A) = P M k=1 p M k M k (A), for A2 where the nodesf M k g M k=1 and the weightsfp M k g M k=1 R + with P M k=1 p M k = 1. In the case of the mixed eects model, the approximating likelihood to be maximized takes the form 105 L N;M (P M ;V ) = m Y i=1 Z L N i (U();V )dP M () = m Y i=1 1 X d=0 Z 1 0 ni Y j=0 '(V j;i y N j;i ())dP M () = m Y i=1 1 X d=0 Z 1 0 ni Y j=0 '(V j;i y N j;i ())d M X k=1 p M k M k () = m Y i=1 M X k=1 p M k ni Y j=0 '(V j;i y N j;i ( M k )): In the case of the naive pooled data estimator, we have that v N;M (t j ;P M ) =E[y(t j ;x N 0 ; u;U)jP M ] = Z y N j (U())dP M () = 1 X d=0 Z 1 0 y N j (U())d M X k=1 p M k M k () = M X k=1 p M k y N j ( M k ); (12.0.4) in which case J N;M (P M ;V ) = n X j=0 (V j v N;M (t j ;P M )) 2 = n X j=0 (V j M X k=1 p M k y N j ( M k )) 2 : (12.0.5) The optimization can be with respect just thefp M k g M k=1 with thef M k g M k=1 chosen a-priori and xed, or the optimization can be with respect to both the fp M k g M k=1 and thef M k g M k=1 . In the case where the statistical model is the naive pooled data estimator and the optimization is over only thefp M k g M k=1 , the op- 106 timization problem given in (12.0.4) and (12.0.5) becomes one of minimizing a quadratic form over a subset of M dimensional Euclidean space, ^ P N;M = arg min P M 2P() J N;M (P M ;V ) = arg min p M 2R + M; P M k=1 p M k =1 f(p M ) T Hp M +L T p M g where H k;l = n X j=0 y N j ( M k )y N j ( M l );k;l = 1; 2;:::;M; and L k =2 n X j=0 V j y N j ( M k );k = 1; 2;:::;M: With regard to verifying the hypotheses of the convergence theorem for the Prohorov framework in the previous chapter, in order to verify condition (2), it need only be argued that y N j ()!y j ();N!1; uniformly in for 2 . Alternatively, we could use the Prohorov framework to simply estimate the distribution ofU directly. In this case, we assume thatR(U = [0; ^ u], where ^ u is chosen suciently large. We would then let P M (A) = P M k=1 p M k u M k (A), for A2 R(U) where the nodesfu M k g M k=1 R(U) and the weightsfp M k g M k=1 R + with P M k=1 p M k = 1. In this case, (11.3.2) - (11.3.7) become (12.0.6) - (12.0.11) L N;M (P M ;V ) = m Y i=1 M X k=1 p M k ni Y j=0 '(V j;i y N j;i (u M k )); (12.0.6) 107 v N;M (t j ;P M ) = M X k=1 p M k y N j (u M k ); (12.0.7) in which case J N;M (P M ;V ) = n X j=0 (V j M X k=1 p M k y N j (u M k )) 2 ; (12.0.8) and as before ^ P N;M = arg min P M 2P(R(U)) J N;M (P M ;V ) = arg min p M 2R + M; P M k=1 p M k =1 f(p M ) T Hp M +L T p M g; (12.0.9) but now with H k;l = n X j=0 y N j (u M k )y N j (u M l );k;l = 1; 2;:::;M; (12.0.10) and L k =2 n X j=0 V j y N j (u M k );k = 1; 2;:::;M: (12.0.11) 108 Chapter 13 Quantitatively Characterizing and Estimating Drink- ing Patterns and Behavior 13.1 Interpretation of Drinking Parameters In this chapter, we describe how we implemented the scheme presented in Chapter 12 for the estimation of drinking patterns from TAC data based on our random discrete time semi-linear dynamical system modeling alcohol metabolism and transport in the human body that was developed earlier in the thesis. In this chapter we also present the results of some of our numeri- cal experiments and studies. There are two basic components to the numerical experiments we discuss: (i) the simulation of data using a stochastic model de- scribed in Section 13.2 and (ii) using the simulated data or experimental data to estimate the drinking parameters that were used to generate the simulated data or underlying drinking parameters for given subjects. In this way, we be- lieve that we will have demonstrated that our approach yields a reasonable way to quantitatively characterize drinking patterns and behavior. We note, how- ever, that we are only proposing one possible way to do this and that there are of course other ways that our general approach and schemes could be used to quantitatively describe drinking behavior. In the series of experiments to follow, we characterize drinking behavior using two random drinking parameters, d and T. We assume that the number of drinks ingested by the subject during a drinking episode is a random parameter d. Then we assume that duration of a drinking episode is another random 109 parameterT. The numerical study we present in this chapter is to estimate the distribution of the two random drinking parametersd andT. In Section 13.3, we assume thatd andT have parametric density functions, and we can estimate the distributions of d and T by estimating the density parameters. We assume thatd is Poisson distributed with parameter . Ac- cording to the property of Poisson distribution, also represents the expected number of drinks consumed by the subject during a drinking episode. Then we assume that the duration of a drinking episode T has a gamma distribution withT Gamma(;). According to the properties of the Gamma distribu- tion, we know that the expected duration of a drinking episode is then given by =. In Section 13.4, we assume that the distributions of d and T are non- parametric, and we will use the Prohorov framework we developed in Chapter 12 to estimate the distributions ofd andT. 13.2 Data Simulation In this section, we discuss the simulator we used to generate TAC data based on a particular drinking diary as input. In order to keep the simulated data independent of the estimator, we will use a somewhat dierent model for alcohol metabolism and transport than the one upon which our estimation scheme is based. This section has two parts: rst, we provide a brief overview of the simulation model and point out the dierences between this model and the model our scheme is based on. We will also discuss the dierences between the particular numerical methods used in the two models. In the second part, we describe how the simulator actually works and present a few samples of BAC/BrAC and TAC data generated using the simulator. 110 13.2.1 The Stochastic Model There are two essential components that form the underlying basis for the simulator model. The rst is a linear input/output model for TAC, y, as a function of BrAC, w, in the form of a forward convolution y(t) = Z t 0 k(q;t)w()d; (13.2.1) where the unknown convolution kernel,k, is parametrized by a nite dimensional vector of parameters, q2Q, Q a set of feasible parameters. One way to obtain the above input/output model is to consider a one di- mensional diusion equation to relate BrAC and TAC by assuming ethanol molecules from the blood diuse through the epidermal layer of the skin. Thus, the following partial dierential equation can be used to derive the convolution form as given in (13.2.1): @' @t (t;x) =q 1 @ 2 ' @x 2 (t;x); 0x 1;t> 0 q 1 @' @x (t; 0) ='(t; 0);t> 0 q 1 @' @x (t; 1) =q 2 w(t);t> 0 '(0;x) = 0; 0x 1 y(t) ='(t; 0): (13.2.2) To obtain y in terms of a convolution with w from (13.2.2), the above PDE is reformulated in weak form and then operator semi-group theory is used to obtain the solution to the weak form of the of initial boundary value problem given in (13.2.2). This is basically very similar to the process we used to solve our semilinear hybrid PDE-ODE system as was explained in earlier chapters. Consequently, we omit this part of the discussion here. For additional details, 111 one can see [13, 27] as we will proceed directly to the conclusion they drew. To see how the convolution kernel k(q;) can be parameterized by q2 Q via the PDE model given in (13.2.2), we rely on the treatments in [29, 31]. Let Q be a compact subset of R 2 endowed with the Euclidean metric, let H = L 2 (0; 1) together with the standard inner producth 1 ; 2 i = R 1 0 1 2 d, and norm denoted byjj, and letV be the Sobolev spaceH 1 (0; 1) together with its standard inner producthh 1 ; 2 ii = R 1 0 1 2 d + R 1 0 0 1 0 2 d and norm denoted byjjjj. The usual dense and continuous embeddingsV ,!H ,!V hold, where V =H 1 (0; 1) denotes the dual space of V . Dene the forms and functionals a(q;;) :VV !R, b(q;) :V !R, and c(q;) :V !R by a(q; 1 ; 2 ) = 1 (0) 2 (0) +q 1 Z 1 0 0 1 0 2 d; (13.2.3) for 1 ; 2 2V ,hb(q); i V ;V =q 2 (1), andhc(q); i V ;V = (0), for 2V . It follows thatb(q) =q 2 (1)2V andc(q) =2V , where denotes the Dirac delta distribution, or unit impulse at zero. It is not dicult to argue that the forma(;;) as given in (13.2.3) above is bounded, coercive and continuous with repect toq2Q (see, for example, [38], also, see [12] for a more abstract, detailed and rigorous description of how to deal with input signals on the boundary of the domain). Then, forq2Q, the bilinear forma(q;;), denes a bounded linear operator A(q)2L(V;V ) byhA(q) 1 ; 2 i V ;V =a(q; 1 ; 2 ), for 1 ; 2 2V . Then, if we letH denote any of the spacesV;H orV , we can consider the linear operator A(q) to be the unbounded operator, A(q) :D q H!H whereD q =V in the caseH =V , and D q =f 2V :A(q) 2Hg in the caseH =H orH =V . It can be shown (see, for example, [3, 4, 38]) that A(q) is a closed, densely dened unbounded linear operator on H and it is the innitesimal generator of an analytic semigroup of bounded linear operators,fe A(q)t : t 0g onH. 112 Next dene the bounded linear operators B(q) :R!V and C(q) :V !R by hB(q)u; i V ;V =hb(q); i V ;V u and C(q) =hc(q); i V ;V , respectively, for 2V and u2R. y(t) = Z t 0 Ce A(q)(ts) B(q)w(s)ds; whereA(q)2L (V;V ) is a bounded linear operator dened by<A(q) 1 ; 2 >= a(q; 1 ; 2 ). It follows that k(q;t) =Ce A(q)t B(q). The second component of the simulator is the nonlinear pharmacokinetic model for alcohol absorption based on Michaelis-Menten enzyme kinetics _ w(t) =bu(t) Kw(t) M +w(t) ;t> 0;w(0) = 0; (13.2.4) wheret is time in hours(h),w(t) is the subject's BAC or BrAC in units of mg/dl, u(t) is the number of standard drinks (sds) ingested by the subject at time t. The simulator uses the above model, (13.2.4) and (13.2.1) to generate simu- lated drinking episodes as follows. First, a drink dairy is formulated. Then the number of standard drinks consumed at timet is substituted foru(t) in (13.2.4). This ODE is integrated using a standard 4th order Runge-Kutta method to ob- tain BAC or BrAC, w(t). Then the BAC or BrAC, w(t) is substituted into (13.2.1) to obtain TAC,y(t). Note that the primary dierence between the sim- ulator and our model is that the simulator consists of two distinct equations, whereas our model is a hybrid PDE-ODE system. Therefore, while we solve the PDE-ODE system as one system in discrete time, the simulator consists of a two step process in continuous time. The parameters in the simulator were obtained by tting the two step model to actual drink diary.BrAC/TAC data (see [13]). 113 13.2.2 The Simulator and Sampled Data As explained in the last section, there is essentially two steps involved in using the simulator. The simulator was implemented in MATLAB and Simulink. One drinking episode of sampled data generated by the simulator consists of 1) a drinking diary, 2) BAC or BrAC and TAC data over a period of time with a predened sampling interval of hours. And we also designated the duration of the data recording to beT R hours. So the number of sampled points for each drinking episode isn i = T R . By repeating the data simulation processm times, we will end up with m drinking episodes where for the i-th drinking episode we have drinking diary data as a 1 2 vector, and BAC or BrAC and TAC data each as a 1n i vector. We will show in the following (1) how the drinking diary is generated, and (2) how BrAC and TAC are computed using the model described in the previous section. To simulate a drinking diary so that it re ects the way that an actual drink- ing diary is generated, we assume that the subject consumes d i number of standard drinks in the i-th drinking episode which lasts T i hours. We make the assumption (see Chapter 12) that d i and T i are random and sampled from d and T, respectively. We used the built-in funtion POISSRND with given parameter in MATLAB to generate d i by simply calling "poissrnd()" in MATLAB. Similarly, we used the built-in funtion GAMRND in Matlab with given parameters and to generate T i by simply calling "gamrnd(;)" in MATLAB. So far, (d i ;T i ) is used to generate the drinking diary for the i-th drinking episode. The inputu i in the nonlinear model (13.2.4) will then be a 1n i vector with the rst j Ti T R n i k elements beingd i =T i and the rest of the elements being zero. In order to get BrAC and TAC, we assume that there exists either a set of 114 population parameters q 1 ;q 2 for the PDE model (13.2.2) and ;K;M for the ODE model (13.2.4), or parameters that have been t to an individual subject. In th examples that follow, we used the set of parameter values q 1 = 0:99;q 2 = 0:58;b = 0:76;K = 0:90; and M = 2:09. These parameter values were obtained from our numerical studies in Section 10.6. Recall that this set of parameters were estimated by tting drinking episodes 1-7 from the dataset described in Section 10.3. Using a standard 4th order Runge-Kutta method applied to (13.2.4), we can obtain a discrete approximation, w i , to w which is BAC or BrAC data for the i-th drinking episode. Then the BAC or BrAC, w i , is substituted into (13.2.1) to obtain TAC, y i . Note that the scheme we used to solve (13.2.1) numerically is omitted here, since we are using the implementation in [27, 13]. In the fol- lowing plots, we show two illustrative drinking episodes with dierent values for the parameters. Both of the drinking episodes were simulated under the as- sumption that BrAC and TAC measurements were collected over a period of 15 hours. The drinking episode below was generated from a random drinking diary generated using Poisson distribution parameter = 5, and Gamma distribution parameters = 1 and = 1. The next two gures are the BrAC and the TAC vs time. We repeated the above simulation process 20 times yielding 20 drinking episodes, each with a drinking diary, and BAC or BrAC and TAC curves. We will use this set of sampled data in the next section to show how the drinking parameter estimation problem will be done. 115 Figure 13.2.1: Simulated BrAC with density parameters = 5; = 1; = 1 Figure 13.2.2: Simulated TAC with density parameters = 5; = 1; = 1 13.3 TheParameterEstimationSchemeforPara- metricallyEstimatingtheProbabilityDen- sity Function First, we want to be clear about the parameter estimation problem we are trying to solve. Suppose that the training data we have is TAC measurements fV j;i g, j = 1; ;n i ;i = 1; ;m. In other words, there are m drinking 116 episodes, and there are n i TAC measurements in the i-th drinking episode. As was discussed in Chapter 12, we have two dierent methods to estimate drinking parameters: (1) We can parametrically estimate the joint density func- tion of the drinking parametersd andT, or (2) We can use the Prohorov metric framework to estimating the distribution of the drinking parametersd andT directly in the form of probability measure. We only discuss the rst method in this section. The second method will be discussed in the next section. In addi- tion, for each estimation method, we consider two dierent statistical models as was discussed in Chapter 12, the mixed eect model and the naive pooled data model. We have carried out numerical experiments for each one. For the mixed eects model, we have the following optimization problem. For simplicity, we formulate the parameter estimation problem as a maximum likelihood estimation problem and, by taking negative logarithm of the likeli- hood (9.2.1), we obtain the following minimization problem ^ N = arg min 2 logL N (f ;V ) = arg min 2 m X i=1 log 1 X d=0 Z 1 0 ni Y j=0 '(V j;i y N j;i )f (d;T )dT; (13.3.1) where y N j;i is given by ' N j+1;i =S N (;q)' j;i + ^ B N (q)F (i;' N j;i ;u j;i ;d=T;q) ' N 0;i = 0 y N j;i =C(q)' N j;i ;j = 1; 2; ;n i ;i = 1; 2; ;m; (13.3.2) For the naive pooled data model, we formulate the parameter estimation 117 problem as ^ N = arg min 2 J N (f ;V ) = arg min 2 n X j=0 (V j v N (t j ;f )) 2 : (13.3.3) where v(t j ;) =E[y j (U)jf ] = 1 X d=0 Z 1 0 y j (U)f (d;T )dT; represents the mean behavior across all drinking episodes at timet j ,j = 0;:::;n and y N j is given by (13.3.2). To solve the above optimization problems, we used the MATLAB built-in routine PATTERNSEARCH in the MATLAB Stochastic Optimization Toolbox. As a derivative-free search, PATTERNSEARCH needs two inputs: a function handle and an initial value. The function handle computes the value of the objective function. The initial value is a real vector specifying an initial guess for the parameters for the pattern search algorithm. The basic idea of a pattern search algorithm is to nd a sequence of points that approach the optimal so- lution. The scheme starts with the initial guess x 0 , and then adds the pattern vectors to the initial point with a given mesh size to get mesh points. The algorithm evaluates the values of the objective funtion at the mesh points, and nds the location of the minimum among all the function values, denoted by x 1 . If the value of the function is smaller than the value at x 0 , then this is a successful poll, and the algorithm sets the next point in the sequence equal to x 1 and then enlarges the mesh size by 2. If the function value at x 1 is not smaller than the value atx 0 , then this is an unsuccessful poll and the algorithm then shrinks the mesh size by 0.5. The algorithm stops when the mesh size is less than a specied minimum mesh tolerance (the default in MATLAB is 1e-6). 118 Note that the actual working of a pattern search algorithm is not our primary concern here, we are just just using it as a global search tool. The reason why we chose a pattern search algorithm as our optimization tool is that a pattern search algorithm does not require the gradient of the objective functions given in (13.3.1) and (13.3.3) and because it is relatively ecient when compared to other Matlab built-in optimization routines. The routine PATTERNSEARCH requires that we supply it with a function that evaluates the objective functions as given in (13.3.1) and (13.3.3). This raises two questions: (1) how should y N j;i in system (13.3.2) be computed using a specic nite dimensional approximation, and then (2) how exactly should we evaluate the objective functions in (13.3.1) and (13.3.3). The rst question was answered in Section 10.2. In order to evaluate L N (f ;V ) in (13.3.1), we apply a Monte Carlo inte- gration technique. The basic idea of a Monte Carlo integration technique is to generate a number of samplesf(d;T )g L k=1 according to a joint density funtion (12.0.3) with density parameters , and . So the k-th drinking sample lasts T i hours and d i standard drinks are ingested during the k-th drinking episode. The length of the sampling interval is hours; consequently u k =d k =T k stan- dard drinks per hour. If we let k = (T k ;d k ) and u j;k = , for all k;j, then the amount of alcohol ingested in each sampling interval in units of standard drinks is given byu k u j;k =u( k )u j;k = (d k =T k )u j;k = (d k =T k ) standard drinks. Substituting u k u j;k into 13.3.2, we obtain y N j;k . We then compute ni Y j=0 '(V j;i y N j;k ) 119 for each (d k ;T k ). A crude way of calculating the integral is L N (f ;V ) = m X i=1 log[ 1 L L X k=1 ni Y j=0 '(V j;i y N j;k )]: (13.3.4) In order to evaluate J N (f ;V ) in (13.3.3), we require v i (t j ;). Using the same Monte Carlo scheme to generatey N j;k , we apply a Monte Carlo integration method and obtain v(t j ;) = 1 L L X k=1 y N j;k : (13.3.5) 13.4 The Parameter Estimation Scheme for Es- timatingtheProbabilityMeasureBasedon the Prohorov Metric Framework To use the Prohorov metric framework approach for estimating the distribu- tion of random parameters in dynamical systems, we assume that the duration of drinking,T, is known, so the only random drinking parameter isd. We then used the Prohorov framework to estimate the distribution ofd. Let D =f0; 1; 2; 3;:::; ^ dg, where ^ d is suciently large and xed (recall, we wantD to be compact). Then as we noted in Chapter 12, the space of probabil- ity measuresP(D) can be discretized by considering measures that are linear combinations of Dirac measures of the form P M (A) = P M k=1 p M k M k (A), for A2 D where the nodesfd M k g M k=1 D and the weightsfp M k g M k=1 R + with P M k=1 p M k = 1. With the number of nodes M and the choice of nodesfd M k g M k=1 specied, we seek the optimal weights P M =fp M k g M k=1 . In the case of the mixed eects model, the approximating negative log like- 120 lihood to be minimized takes the form ^ P M = arg min P M logL N;M (P M ;V ) = arg min P M m X i=1 log M X k=1 p M k ni Y j=0 '(V j;i y N j;i (d M k )): (13.4.1) In the case of naive pooled model, the optimization problem is ^ P M = arg min P M J N;M (P M ;V ) = arg min P M n X j=0 (V j M X k=1 p M k y N j (d M k )) 2 : (13.4.2) To solve the above optimization problems, once again we use the MATLAB built-in routine PATTERNSEARCH in the MATLAB Stochastic Optimization Toolbox that was described in the previous section. Evaluate the objective funtion in this case is really quite straight forward. Fix the value for M and dene the values of nodesfd M k g M k=1 . For any given P M =fp M k g M k=1 , we can evaluate the objective functions in (13.4.1) and (13.4.2) directly since the they N j;i (d M k )'s can be computed using the same scheme as the one described in Section 10.2. 13.5 Numerical Results 13.5.1 Numerical Results for Parametrically Estimating the Probability Density Function In the numerical experiments described in this section, we set the length of the sampling interval to be = 1=20 hour in the system (12.0.2). The number of spatial statesN in the nite dimensional approximating systems (13.3.2) was set to 4. The number of Monte Carlo samples L in the approximation formulas 121 (13.3.4) and (13.3.5) was taken to be 100. In the two numerical experiments that follow, We chose the initial guess 0 = 1, and we imposed the compactness constraint be requiring that 2 [0:1; 10]. The optimization problem was solved using the routine PATTERNSEARCH in the Matlab Stochastic Optimization Toolbox. 13.5.2 Numerical Experiment 1 Data: We assumed that the duration of drinking during the episode is xed atT = 5. We chose the density parameter to be = 5, and we generatedfd i g 20 i=1 using the MATLAB routine POISSRND to obtain 20 Poisson distributed random drinks per episode. Then using the simulator described in Section 13.2, 20 drinking episodes of TAC datafV j;i g, i = 1; ; 20, j = 1; ; 300 were generated with the duration of measurement being 15 hours. Results: Based on the data described above with T = 5 and the naive pooled data statistical model as described in Section 13.3, the optimal estimate of the Poisson parameter that was used to generate the data, the solution to the optimization problem (13.3.3) was = 5:0137. 13.5.3 Numerical Experiment 2 Data: Same as in Numerical Experiment 1. Results: Based on the data described above with T = 5 and the naive pooled data statistical model as described in Section 13.3, the optimal estimate of the Poisson parameter that was used to generate the data, the solution to the optimization 122 problem (13.3.3) was = 4:9980. In the next four numerical experiments, we chose the initial guesses = 1, = 1, and = 1. We imposed the compactness constraints be requiring that 2 [0:1; 10], 2 [0:1; 10] and 2 [0:1; 2], when the optimization problem was solved once again using the routine PATTERNSEARCH in the Matlab Stochastic Optimization Toolbox. 13.5.4 Numerical Experiment 3 Data: When we generated the data, we chose the density parameters to be = 5, = 5, and = 1. We generatedfd i ;T i g 20 i=1 using the MATLAB routines POISSRND and GAMRND respectively to generate 20 Poisson samples for the number of drinks per episode and 20 Gamma distributed samples for the duration of drinking during the episode. Then using the simulator described in earlier in Section 13.2, 20 drinking episodes with TAC datafV j;i g,i = 1; ; 20, j = 1; ; 300 were generated with the duration of measurement being taken to be 15 hours. Results: Based on the data described above and the naive pooled data statistical model as described in Section 13.3, the optimal estimates of the Poisson parameter and the Gamma parameters that were used to generate the data, the solution to the optimization problem (13.3.3) was = 5:75 , = 6:0, and = 1:0. 13.5.5 Numerical Experiment 4 Data: Same as in Numerical Experiment 3. Results: 123 Based on the data described above and the mixed eects statistical model as described in Section 13.3, the optimal estimates of the Poisson parameter and the Gamma parameters that were used to generate the data, the solution to the optimization problem (13.3.1) was = 4:03 , = 5:0, and = 1:0. 13.5.6 Numerical Experiment 5 Data: Experimental data described in Section 10.3. Results: Based on the data described above and the naive pooled data statistical model as described in Section 13.3, the optimal estimates of the Poisson parameter and the Gamma parameters that were used to generate the data, the solution to the optimization problem (13.3.3) was = 1:63 , = 1:48, and = 0:50. 13.5.7 Numerical Experiment 6 Data: Experimental data described in Section 10.3. Results: Based on the data described above and the mixed eects statistical model as described in Section 13.3, the optimal estimates of the Poisson parameter and the Gamma parameters that were used to generate the data, the solution to the optimization problem (13.3.1) was = 1:25 , = 1:0, and = 0:88. 124 13.6 Numerical Results for Non-Parametrically Estimating the Drinking Parameters using Prohorov Metric Framework In the numerical experiments described in this section, we set the length of the sampling interval to be = 1=20 hour in the system (12.0.2). The number of spatial statesN in the nite dimensional approximating systems (13.3.2) was set to 4. In the two numerical experiments that follow, the optimization problem was solved using the routine FMINCON in the Matlab Stochastic Optimization Toolbox. In numerical experiment 7, we used a Poisson distribution to generate the number of drinks in each episode, so d is a discrete random variable. In nu- merical experiment 8, we used a Normal distribution to generate the number of drinks in each episode, so in this experimentd is a continuous random variable. 13.6.1 Numerical Experiment 7 Data: We assumed that the duration of drinking during the episode is xed atT = 5. We assumed that the number of drinks d follows a Poisson distribution with parameter = 5, and we generated 1000 samples,fd i g 1000 i=1 , using the MATLAB routine POISSRND to simulate 1000 drinking episodes, the ith with d i drinks per episode. Then using the simulator described in Section 13.2, 1000 drinking episodes of TAC datafV j;i g, i = 1; ; 1000, j = 1; ; 300 were generated with the duration of measurement being 15 hours. Results: 125 Based on the data described above withT = 5 and the mixed eects statistical model as described in Section 13.4, we sought to estimate the weightsfp M k g M k=1 with P M k=1 p M k = 1 inP M (A) = P M k=1 p M k M k (A), forA2 where the nodes f M k g M k=1 are dened as [0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10] and M = 10. we chose the initial guess p k = 1=11;k = 1; ; 11, and we imposed the compactness constraint by requiring that 0 p k 1;k = 1; ; 11 and P 11 k=1 p k = 1. The optimal estimate of thep k s obtained as the solution to the optimization problem (13.3.1) was p = [:011;:039;:071;:125;:183;:187;:131;:119;:061;:042;:032]: (13.6.1) The bar graph on the left in Figure 13.6.1 shows our approximation to the distribution ofd as given in (13.6.1), and the bar graph on the right in Figure 13.6.1 shows the actual distribution ofd used to simulate the data. Notice that the approximating distribution is close to the actual distribution, which is the Poisson distribution with = 5. Figure 13.6.1: Left:Approximation of Distribution ofd using Prohorov Metric; Right: Actual Distribution ofd used to Generate Data 126 13.6.2 Numerical Experiment 8 Data: We assumed that the duration of drinking during the episode is xed at T = 5. We assumed that the number of drinks d follows a Normal distribution with = 5 and = 1, and we generated 1000 samples,fd i g 100000 i=1 , using the MATLAB routine NORMRND to obtain 10000 drinking episodes with normally distributed number of drinks per episode. Then using the simulator described in Section 13.2, 100000 drinking episodes of TAC datafV j;i g,i = 1; ; 100000, j = 1; ; 300 were generated with the duration of measurement being 15 hours. Results: Based on the data described above withT = 5 and the mixed eects statisti- cal model as described in Section 13.4, we sought to approximate the weights fp M k g M k=1 R + with P M k=1 p M k = 1 inP M (A) = P M k=1 p M k M k (A), forA2 where the nodesf M k g M k=1 are dened as [2; 2:5; 3; 3:5; 4; 4:5; 5; 5:5; 6; 6:5; 7; 7:5; 8] andM = 13. we chose the initial guessp k = 1=13;k = 1; ; 13, and we imposed the compactness constraint by requiring that 0p k 1;k = 1; ; 13 and P 13 k=1 p k = 1. The optimal estimate of the p k s obtained as the solution to the optimization problem (13.3.1) was p = [0:0004;:0048;:0151;:0595;:1131;:1921; :2312;:1854;:1243;:0508;:0199;:0017;:0016]: (13.6.2) The bar graph on the left in Figure 13.6.2 shows our approximation to the distribution ofd as given in (13.6.2), and the bar graph on the right in Figure 13.6.2 shows the actual distribution ofd used to simulate the data. Notice that the approximating distribution is close to the actual distribution, which is the Normal distribution with = 5; = 1. 127 Figure 13.6.2: Left:Approximation of Distribution ofd using Prohorov Metric; Right: Actual Distribution ofd used to Generate Data 128 Chapter 14 Summary, Concluding Remarks and Discussion of Fu- ture Research 14.1 Summary In this thesis we have developed a semi-linear distributed parameter model for the metabolism and transport of alcohol in the human body and its mea- surement by a transdermal alcohol biosensor. The model is a hybrid partial and ordinary reaction-diusion equation. We use non-linear Michaelis-Menten dynamics to describe the enzyme catalyzed reaction in the liver where ethanol is metabolized. We use a one dimensional diusion equation to model the trans- port of ethanol through the skin. The input to the model is standard drinks consumed and the output is both blood and transdermal alcohol concentration. For the purpose of analysis we generalize the model to an abstract semilinear parabolic system set in a Gelfand triple of Hilbert spaces. We rely on results from linear semigroup theory to argue well-posedness and convergence of nite dimensional approximations. We investigated the estimation of the parameters that appear in the model using a nonlinear least squares approach and based on transdermal measurements. We then considered the estimation of random parameters appearing in the model and looked at two dierent approaches: one based on the Prohorov metric on the space of probability measures and one based the estimation of the probability density function for the random pa- rameters considered as random variables. We then showed how our approach could be used to quantitatively characterize and estimate drinking patterns and 129 behavior based on measurements of transdermal alcohol. 14.2 Concluding Remarks In this thesis, we established a method to quantitatively characterize drink- ing patterns and behaviour from TAC data collected on wearable transdermal alcohol biosensors. This is a non-intrusive, passive and somewhat reliable way to characterize drinking patterns in related clinical studies. If a transdermal alcohol biosensor could potentially become a part of a wearable vital sign mon- itoring device (for example, like a Fitbit), then our method could be used to monitor the wearer's drinking behaviour and provide data and insight that could be used in clinical drinking pattern studies across a cohort or population. 14.3 Avenues of Possible Future Research Several possible future research problems remain to be considered. First, the semi-linear distributed parameter model for the metabolism and transport of alcohol in this thesis really only includes tow organs, the liver and the skin. In fact, the transport and metabolism of alcohol in the human body involves and more signicantly, aects. numerous other organ systems including the stomach, the gut, and especially, the brain. It would be interesting to include these other sites in the body as compartments in our model and to see 1) if any associated parameters could be estimated based on measurements of transdermal alcohol concentration, and 2) what if anything could be learned about how alcohol aects these organ systems from the values of the t parameters. mathematical questions that remain include several involving the very im- portant notion of the identiability of the parameters. Unfortunately for a system of equations as complex as the one we have studied here, such an inves- 130 tigation would likely be challenging if not impossible. Additionally, for the com- putational scheme based on the Prohorov metric framework, it is not clear how to choose the number of nodesM and the values of the nodesf k g M k=1 from the dense subset of . Theoretically, as long as the nodesf M k g M k=1 cover a reason- able subset of the space , this should be adequate. But a more in-depth study is probably warranted. In addition, it would be interesting to look at schemes in which both the nodes,f M k g M k=1 and the weightsfp M k g M k=1 are optimized. This of course would lead to a signicantly more complex and computationally challenging optimization problem, especially in the case of the mixed eects statistical model. Finally, we are interested in considering Bayesian schemes for estimating both the model parameters and the drinking parameters (see, for example, [8, 9, 10, 14, 32, 33]). 131 Bibliography [1] Jones AW. Evidence-based survey of the elimination rates of ethanol from blood with applications in forensic casework. Forensic Sci. Int., 200(1{3):1{ 20, 2010. [2] H. T. Banks, K. B. Flores, I. G. Rosen, E. M. Rutter, Melike Sirlanci, and Clayton Thompson. The prohorov metric framework and aggregate data inverse problems for random pdes. Communications in Applied Analysis, 22(3):415{446, 2018. URL: https://acadsol.eu/en/articles/22/3/6. pdf, doi:10.12732/caa.v22i3.6. [3] H. T. Banks and K. Ito. A unied framework for approximation in inverse problems for distributed parameter systems. Control Theory Advanced Technology, 4(1):73{90, 1988. URL: https://apps.dtic.mil/dtic/tr/ fulltext/u2/a193780.pdf. [4] H. T. Banks and Karl Kunisch. Estimation Techniques for Distributed Parameter Systems. Springer Science & Business Media, 2012. URL: https://www.springer.com/us/book/9780817634339. [5] H. T. Banks and W Clayton Thompson. Least squares estimation of proba- bility measures in the prohorov metric framework. Technical report, DTIC Document, 2012. URL: https://www.researchgate.net/publication/ 268353806_Least_Squares_Estimation_of_Probability_Measures_ in_the_Prohorov_Metric_Framework. [6] H.T. Banks and K. Ito. Approximation in lqr problems for innite dimen- sional systems with unbounded input operators. J. Mathematical Systems Estimation and Control, 7(1):1{34, 1997. 132 [7] Bankole Johnson Boris Kovatchev, Marc Breton. In silico models of alcohol dependence and treatment. Front Psychiatry, 3, 2012. [8] Tan Bui-Thanh, Carsten Burstedde, Omar Ghattas, James Martin, Georg Stadler, and Lucas C. Wilcox. Extreme-scale UQ for Bayesian in- verse problems governed by PDEs. In SC12, November 10-16. IEEE, 2012. URL: https://ieeexplore.ieee.org/document/6468442, doi: 10.1109/SC.2012.56. [9] Tan Bui-Thanh, Omar Ghattas, James Martin, and Georg Stadler. A computational framework for innite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion. SIAM J. Sci. Stat. Comp., 35(6):A2494{A2523, 2010. URL: https://epubs.siam.org/doi/abs/10.1137/12089586X? journalCode=sjoce3, doi:https://doi.org/10.1137/12089586X. [10] Daniela Calvetti, Jario P. Kaipio, and Erkki Somersalo. Inverse problems in the Bayesian framework. Inverse Problems, 30:1{4, 2014. URL: iopscience.iop.org/article/10.1088/0266-5611/30/11/ 110301/pdf, doi:10.1088/0266-5611/30/11/110301. [11] George Casella and Roger L Berger. Statistical Inference, volume 2. Duxbury, Pacic Grove, CA, 2002. URL: https://books.google.com/ books/about/Statistical_Inference.html?id=0x_vAAAAMAAJ. [12] R. F. Curtain and D. Salamon. Finite dimensional compensators and in- nite dimensional systems with unbounded input operators. SIAM J. Con- trol and Optimization, 24(4):797{816, 1986. [13] Zheng Dai, I Gary Rosen, Chunming Wang, Nancy Barnett, and Susan E Luczak. Using drinking data and pharmacokinetic modeling to calibrate 133 transport model and blind deconvolution based data analysis software for transdermal alcohol biosensors. Mathematical Biosciences and Engineering: MBE, 13(5):911, 2016. URL: http://www.aimsciences.org/journals/ displayArticlesnew.jsp?paperID=12739. [14] Masoumeh Dashti and Andrew M. Stuart. The Bayesian approach to in- verse problems. In R. Ghanem et al., editor, Handbook of Uncertainty Quantication, pages 311{428. Springer International Publishing Switzer- land, 2017. URL: https://www.springer.com/us/book/9783319123844. [15] M. Davidian and D. Giltinan. Nonlinear Models for Repeated Measure- ment Data. Chapman and Hall, New York, 1995. URL: https://www. crcpress.com/Nonlinear-Models-for-Repeated-Measurement-Data/ Davidian-Giltinan/p/book/9780412983412. [16] M. Davidian and D. M. Giltinan. Nonlinear models for repeated mea- surement data: An overview and update. Journal of Agricultural, Bi- ological and Environmental Statistics, 8:387{419, 2003. URL: https: //link.springer.com/article/10.1198/1085711032697. [17] E. Demidenko. Mixred Models, Theory and Applications, Sec- ond Edition. John Wiley and Sons, Hoboken, 2013. URL: https://www.wiley.com/en-us/Mixed+Models%3A+Theory+and+ Applications+with+R%2C+2nd+Edition-p-9781118091579. [18] Miguel A Dumett, I Gary Rosen, J Sabat, A Shaman, L Tempelman, C Wang, and R. M. Swift. Deconvolving an estimate of breath measured blood alcohol concentration from biosensor collected transdermal ethanol data. Applied Mathematics and Computation, 196(2):724{743, 2008. URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2597868/, doi:10. 1016/j.amc.2007.07.026. 134 [19] Claude Jerey Gittelson, Roman Andreev, and Christoph Schwab. Op- timality of adaptive galerkin methods for random parabolic partial dif- ferential equations. J. Computational Applied Mathematics, 263:189{ 201, 2014. URL:http://dx.doi.org/10.1016/j.cam.2013.12.031,doi: 10.1016/j.cam.2013.12.031. [20] Tosio Kato. Perturbation Theory for Linear Operators, volume 132. Springer Science & Business Media, 2013. URL: https://www.springer. com/us/book/9783540586616. [21] Dominick A. Labianca. The chemical basis of the breathalyzer: a criti- cal analysis. J. Chem. Educ, 67(3):259{261, 1990. URL: https://pubs. acs.org/doi/abs/10.1021/ed067p259?journalCode=jceda8, doi:DOI: 10.1021/ed067p259. [22] A FJ Levi and I Gary Rosen. A novel formulation of the adjoint method in the optimal design of quantum electronic devices. SIAM Journal on Control and Optimization, 48(5):3191{3223, 2010. URL: https://doi. org/10.1137/070708330, arXiv:https://doi.org/10.1137/070708330, doi:10.1137/070708330. [23] E Nyman and A Palml ov. The elimination of ethyl alcohol in sweat. Acta Physiologica, 74(2):155{159, 1936. URL: https://onlinelibrary. wiley.com/doi/abs/10.1111/j.1748-1716.1936.tb01150.x, doi:https://doi.org/10.1111/j.1748-1716.1936.tb01150.x. [24] Widmark E. M. P. Die theoretischen Grundlagen und die praktische Ver- wendbarkeit der gerichtlich-medizinischen Alkoholbestimmung. Urban and Schwarzenberg, 1932. 135 [25] A. Pazy. Semigroups of Linear Operators and Applications to Partial Dif- ferential Equations. Applied Mathematical Sciences. Springer, 1983. URL: https://books.google.com/books?id=80XYPwAACAAJ. [26] I. G. Rosen. discrete approximation framework for hereditary systems. Technical report, Ph.D. thesis, Brown Univ.,Providence, RI, 1980. [27] I Gary Rosen, Susan E Luczak, and Jordan Weiss. Blind deconvolu- tion for distributed parameter systems with unbounded input and output and determining blood alcohol concentration from transdermal biosensor data. Applied Mathematics and Computation, 231:357{376, 2014. URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3972634/, doi:10. 1016/j.amc.2013.12.099. [28] Christoph Schwab and Claude Jerey Gittelson. Sparse ten- sor discretizations of high-dimensional parametric and stochas- tic pdes. Acta Numerica, 20:291{467, 2011. URL: https: //www.cambridge.org/core/journals/acta-numerica/article/ sparse-tensor-discretizations-of-highdimensional-parametric-and-stochastic-pdes/ A46BD443A2D1176B448132A271057DCC, doi:10.1017/ S0962492911000055. [29] Melike Sirlanci, S. E. Luczak, and I. G. Rosen. Estimation of the dis- tribution of random parameters in discrete time abstract parabolic sys- tems with unbounded input and output: Approximation and convergence. Communications in Applied Analysis, 23(2):287{329, 2019. URL: https: //doi.org/10.12732/caa.v23i2.4, doi:doi:10.12732/caa.v23i2.4. [30] Melike Sirlanci, Susan Luczak, and I Gary Rosen. Approximation and convergence in the estimation of random parameters in linear holo- morphic semigroups generated by regularly dissipative operators. In 136 American Control Conference (ACC), 2017, pages 3171{3176. IEEE, 2017. URL: https://www.researchgate.net/publication/318333926_ Approximation_and_convergence_in_the_estimation_of_random_ parameters_in_linear_holomorphic_semigroups_generated_by_ regularly_dissipative_operators, doi:10.23919/ACC.2017.7963435. [31] Melike Sirlanci, Susan E. Luczak, Catharine E. Fairbairn, Dayheon Kang, Ruoxi Pan, Xin Yu, and I Gary Rosen. Estimating the distribution of random parameters in a diusion equation forward model for a transdermal alcohol biosensor. Automatica, 2018. to appear, arXiv:1808.04058. URL: https://arxiv.org/abs/1808.04058. [32] Melike Sirlanci, I Gary Rosen, Susan E. Luczak, Catharine E. Fairbairn, Konrad Bresin, and Dayheon Kang. Deconvolving the input to random abstract parabolic systems; a population model-based approach to estimat- ing blood/breath alcohol concentration from transdermal alcohol biosensor data. Inverse problems, 34(12), 2018. arXiv:1807.05088v1. URL: http: //iopscience.iop.org/article/10.1088/1361-6420/aae791/pdf. [33] A. M. Stuart. Inverse problems: A Bayesian perspec- tive. Acta Numerica, pages 451{559, 2010. URL: https://www.cambridge.org/core/journals/acta-numerica/ article/inverse-problems-a-bayesian-perspective/ 587A3A0D480A1A7C2B1B284BCEDF7E23, doi:https://doi.org/10. 1017/S0962492910000061. [34] Robert M. Swift. Transdermal measurement of alcohol consump- tion. Addiction, 88(8):1037{1039, 1993. URL: http://dx.doi.org/10. 1111/j.1360-0443.1993.tb02122.x, doi:10.1111/j.1360-0443.1993. tb02122.x. 137 [35] Robert M. Swift. Transdermal alcohol measurement for estimation of blood alcohol concentration. Alcoholism: Clinical and Experimental Research, 24(4):422{423, 2000. URL: http://dx.doi.org/10.1111/j.1530-0277. 2000.tb02006.x, doi:10.1111/j.1530-0277.2000.tb02006.x. [36] Robert M. Swift. Direct measurement of alcohol and its metabolites. Addic- tion, 98:73{80, 2003. URL: http://dx.doi.org/10.1046/j.1359-6357. 2003.00605.x, doi:10.1046/j.1359-6357.2003.00605.x. [37] Robert M. Swift, Christopher S Martin, Larry Swette, Anthony LaConti, and Nancy Kackley. Studies on a wearable, electronic, transdermal al- cohol sensor. Alcoholism: Clinical and Experimental Research, 16(4):721{ 725, 1992. URL:https://www.ncbi.nlm.nih.gov/pubmed/1530135,doi: https://doi.org/10.1111/j.1530-0277.1992.tb00668.x. [38] H. Tanabe. Equations of Evolution. Monographs and Studies in Math- ematics. Pitman, 1979. URL: https://books.google.com/books?id= Dn6zAAAAIAAJ. [39] Singh P. Fogler H. S Umulis D. M., G urmen N. M. A physiologically based model for ethanol and acetaldehyde metabolism in human beings. Alcohol, 35(3):1210, 2005. [40] Abdelmegeed MA Banerjee A Morgan TR Yoo S-H Song B-J Yun J-W, Son M-J. Binge alcohol promotes hypoxic liver injury through a cyp2e1{hif-1- dependent apoptosis pathway in mice and humans. Free Radic Biol Med, 77:94{183, 2014. 138

Abstract (if available)

Abstract We developed a deterministic semilinear parabolic system to describe alcohol metabolism and transport in the human body based on transdermal alcohol concentration. The novelty of this model is that we used an ordinary differential equation to describe the elimination process of blood alcohol and a partial differential equation to describe the transport of alcohol molecules through the skin. In this model, we assume that ingested alcohol is immediately absorbed into the blood and the amount of alcohol excreted through the skin is proportional to the amount of alcohol in the blood. We also built a computational scheme to train individual's parameters in the alcohol metabolism and transport model from one or multiple drinking episodes of transdermal alcohol concentration data. With trained parameters, we can evaluate blood alcohol concentration from transdermal alcohol concentration. ❧ Later in this thesis, we constructed a random discrete time dynamical system by assuming input parameters in the deterministic semilinear parabolic system are random. Therefore, the drinking behavior estimation problem is formulated as a problem of identifying random parameters in such system. By establishing a Prohorov metric framework to non-parametrically estimate the distribution of the random parameters and a parametric framework to estimate the distribution of the random parameters, we provided a possible way to quantitatively study people's drinking behaviour from transdermal alcohol concentration collected from non-intrusive wearable biosensor.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Finite dimensional approximation and convergence in the estimation of the distribution of, and input to, random abstract parabolic systems with application to the deconvolution of blood/breath al...

PDF

Prohorov Metric-Based Nonparametric Estimation of the Distribution of Random Parameters in Abstract Parabolic Systems with Application to the Transdermal Transport of Alcohol

PDF

Linear quadratic control, estimation, and tracking for random abstract parabolic systems with application to transdermal alcohol biosensing

PDF

An abstract hyperbolic population model for the transdermal transport of ethanol in humans: estimating the distribution of random parameters and the deconvolution of breath alcohol concentration

PDF

Distributed parameter model based system identification and filtering in the estimation of blood alcohol concentration from transdermal alcohol biosensor data

PDF

M-estimation and non-parametric estimation of a random diffusion equation-based population model for the transdermal transport of ethanol: deconvolution and uncertainty quantification

PDF

Obtaining breath alcohol concentration from transdermal alcohol concentration using Bayesian approaches

PDF

Population modeling and Bayesian estimation for the deconvolution of blood alcohol concentration from transdermal alcohol biosensor data

PDF

Determining blood alcohol concentration from transdermal alcohol data: calibrating a mathematical model using a drinking diary

PDF

Simultaneous parameter estimation and semi-blind deconvolution in infinite-dimensional linear systems with unbounded input and output

PDF

A nonlinear pharmacokinetic model used in calibrating a transdermal alcohol transport concentration biosensor data analysis software

PDF

Physics-informed machine learning techniques for the estimation and uncertainty quantification of breath alcohol concentration from transdermal alcohol biosensor data

PDF

Delta Method confidence bands for parameter-dependent impulse response functions, convolutions, and deconvolutions arising from evolution systems described by…

Asset Metadata

Creator Li, Shuang (author)

Core Title Estimation of random input to semi-linear abstract parabolic systems with application to quantitative description of drinking behavior based on trans-dermal alcohol concentration

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Applied Mathematics

Publication Date 07/28/2019

Defense Date 06/05/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag alcohol biosensor,Michaelis-Menten kinetics,naive pooled analysis,nonlinear mixed effects model,OAI-PMH Harvest,Ode,PDE

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Rosen, Gary (committee chair), Luczak, Susan (committee member), Wang, Chunming (committee member)

Creator Email lishuang@usc.edu,shuang.li90@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-196875

Unique identifier UC11663057

Identifier etd-LiShuang-7651.pdf (filename),usctheses-c89-196875 (legacy record id)

Legacy Identifier etd-LiShuang-7651.pdf

Dmrecord 196875

Document Type Dissertation

Format application/pdf (imt)

Rights Li, Shuang

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA