Loading [MathJax]/jax/output/SVG/jax.js
Research article

Investigation of fractal-fractional HIV infection by evaluating the drug therapy effect in the Atangana-Baleanu sense


  • In this paper, we apply the fractal-fractional derivative in the Atangana-Baleanu sense to a model of the human immunodeficiency virus infection of CD4+ T-cells in the presence of a reverse transcriptase inhibitor, which occurs before the infected cell begins producing the virus. The existence and uniqueness results obtained by applying Banach-type and Leray-Schauder-type fixed-point theorems for the solution of the suggested model are established. Stability analysis in the context of Ulam's stability and its various types are investigated in order to ensure that a close exact solution exists. Additionally, the equilibrium points and their stability are analyzed by using the basic reproduction number. Three numerical algorithms are provided to illustrate the approximate solutions by using the Newton polynomial approach, the Adam-Bashforth method and the predictor-corrector technique, and a comparison between them is presented. Furthermore, we present the results of numerical simulations in the form of graphical figures corresponding to different fractal dimensions and fractional orders between zero and one. We analyze the behavior of the considered model for the provided values of input factors. As a result, the behavior of the system was predicted for various fractal dimensions and fractional orders, which revealed that slight changes in the fractal dimensions and fractional orders had no impact on the function's behavior in general but only occur in the numerical simulations.

    Citation: Jutarat Kongson, Chatthai Thaiprayoon, Apichat Neamvonk, Jehad Alzabut, Weerawat Sudsutad. Investigation of fractal-fractional HIV infection by evaluating the drug therapy effect in the Atangana-Baleanu sense[J]. Mathematical Biosciences and Engineering, 2022, 19(11): 10762-10808. doi: 10.3934/mbe.2022504

    Related Papers:

    [1] Hyung-Chun Lee . Efficient computations for linear feedback control problems for target velocity matching of Navier-Stokes flows via POD and LSTM-ROM. Electronic Research Archive, 2021, 29(3): 2533-2552. doi: 10.3934/era.2020128
    [2] Wei Shi, Xinguang Yang, Xingjie Yan . Determination of the 3D Navier-Stokes equations with damping. Electronic Research Archive, 2022, 30(10): 3872-3886. doi: 10.3934/era.2022197
    [3] Guoliang Ju, Can Chen, Rongliang Chen, Jingzhi Li, Kaitai Li, Shaohui Zhang . Numerical simulation for 3D flow in flow channel of aeroengine turbine fan based on dimension splitting method. Electronic Research Archive, 2020, 28(2): 837-851. doi: 10.3934/era.2020043
    [4] Matthew Gardner, Adam Larios, Leo G. Rebholz, Duygu Vargun, Camille Zerfas . Continuous data assimilation applied to a velocity-vorticity formulation of the 2D Navier-Stokes equations. Electronic Research Archive, 2021, 29(3): 2223-2247. doi: 10.3934/era.2020113
    [5] Jie Zhang, Gaoli Huang, Fan Wu . Energy equality in the isentropic compressible Navier-Stokes-Maxwell equations. Electronic Research Archive, 2023, 31(10): 6412-6424. doi: 10.3934/era.2023324
    [6] Derrick Jones, Xu Zhang . A conforming-nonconforming mixed immersed finite element method for unsteady Stokes equations with moving interfaces. Electronic Research Archive, 2021, 29(5): 3171-3191. doi: 10.3934/era.2021032
    [7] Jiwei Jia, Young-Ju Lee, Yue Feng, Zichan Wang, Zhongshu Zhao . Hybridized weak Galerkin finite element methods for Brinkman equations. Electronic Research Archive, 2021, 29(3): 2489-2516. doi: 10.3934/era.2020126
    [8] Caifeng Liu, Pan Liu . On Liouville-type theorem for the stationary compressible Navier–Stokes equations in R3. Electronic Research Archive, 2024, 32(1): 386-404. doi: 10.3934/era.2024019
    [9] Guochun Wu, Han Wang, Yinghui Zhang . Optimal time-decay rates of the compressible Navier–Stokes–Poisson system in R3. Electronic Research Archive, 2021, 29(6): 3889-3908. doi: 10.3934/era.2021067
    [10] J. Bravo-Olivares, E. Fernández-Cara, E. Notte-Cuello, M.A. Rojas-Medar . Regularity criteria for 3D MHD flows in terms of spectral components. Electronic Research Archive, 2022, 30(9): 3238-3248. doi: 10.3934/era.2022164
  • In this paper, we apply the fractal-fractional derivative in the Atangana-Baleanu sense to a model of the human immunodeficiency virus infection of CD4+ T-cells in the presence of a reverse transcriptase inhibitor, which occurs before the infected cell begins producing the virus. The existence and uniqueness results obtained by applying Banach-type and Leray-Schauder-type fixed-point theorems for the solution of the suggested model are established. Stability analysis in the context of Ulam's stability and its various types are investigated in order to ensure that a close exact solution exists. Additionally, the equilibrium points and their stability are analyzed by using the basic reproduction number. Three numerical algorithms are provided to illustrate the approximate solutions by using the Newton polynomial approach, the Adam-Bashforth method and the predictor-corrector technique, and a comparison between them is presented. Furthermore, we present the results of numerical simulations in the form of graphical figures corresponding to different fractal dimensions and fractional orders between zero and one. We analyze the behavior of the considered model for the provided values of input factors. As a result, the behavior of the system was predicted for various fractal dimensions and fractional orders, which revealed that slight changes in the fractal dimensions and fractional orders had no impact on the function's behavior in general but only occur in the numerical simulations.



    In clinical therapy, medical research, and law enforcement, the breathalyzer, developed by Borkenstein based on a redox reaction and Henry's law [1], is used to measure breath alcohol concentration (BrAC), a surrogate for blood alcohol concentration (BAC). Clinicians and researchers consider it to be reasonably accurate to substitute BrAC for BAC, and in general this continues to be the case across different environmental conditions and across different individuals [1]. Nevertheless, collecting near-continuous BrAC samples accurately (i.e., obtaining a deep lung sample that is not contaminated by any existing alcohol remaining in the mouth) is challenging and often impractical in the field.

    Most of the ethanol, the type of alcohol in alcoholic beverages, that enters the human body, is metabolized by the liver into other products that are then excreted. In addition, a portion of ingested ethanol exits the body directly through exhalation and urination [2] and approximately 1 diffuses through the epidermal layer of the skin in the form of perspiration and sweat. The amount of alcohol excreted in this manner is quantified in the form of transdermal alcohol concentration (TAC). TAC has been shown to be largely positively correlated with BrAC and BAC [3]. However, the precise relationship between TAC and BrAC/BAC is complicated due to a number of confounding physiological, technological and environmental factors including, but not limited to, the skin's epidermal layer thickness, porosity and tortuosity, the process of vasodilation as observed through blood pressure and flow rate, the underlying technology of the particular sensor being used and the ambient temperature and humidity.

    Currently, there are a number of different biosensors based on a variety of analog principles that can measure TAC essentially continuously, passively, unobtrusively, and relatively accurately, and make it available for processing in real time. Some of these devices are already commercially available and more are on the way. Several of these biosensors, like the breathalyzer, rely on relatively standard fuel-cell technology (i.e., converting chemical energy into electricity through redox reactions) to effectively count the number of ethanol molecules that evaporate during perspiration from the epidermal layer of the skin in near-continuous time [4]. Figure 1 shows two of these TAC measuring devices: The WrisTASTM7 developed by Giner, Inc. in Waltham, MA, and the SCRAM CAM (Secure Continuous Remote Alcohol Monitor) developed by Alcohol Monitoring Systems, Inc. (AMS) in Littleton, Colorado.

    Figure 1.  WrisTASTM7 (left) and SCRAM CAM (right) transdermal alcohol biosensors.

    Historically, researchers, clinicians, and the courts have always relied on BrAC or, when available, BAC. Consequently, in order to make TAC biosensors practical and accepted by the alcohol community, reliable and consistent means for converting TAC into equivalent BrAC or BAC must be developed, and this involves challenges that must be dealt with as indicated previously. In the past, our approach to developing a method for converting TAC into BrAC or BAC was based on deterministic methods for estimating parameters in distributed parameter systems, such as those described in [5,6]. Our earlier work along these lines has been reported in, for example, [7,8,9]. In these treatments, a forward model in the form of a one-dimensional diffusion equation based on Fick's law [10] with BrAC as the input and TAC as the output was first calibrated (i.e., fit) using BrAC and TAC data collected from the patient or research subject in the clinic or laboratory during what is known as a controlled alcohol challenge. Then, after the same patient or research subject has worn the TAC sensor in the field for an extended period of time (e.g., days, weeks, or even months), the TAC data is downloaded, and the fit forward model is used to deconvolve the BrAC or BAC input from the observed TAC output.

    In order to eliminate the calibration process, we developed a population model-based approach wherein the parameters in the model were assumed to be random. Then, rather than fitting the values of the parameters themselves, their distributions were estimated based on BrAC and TAC data from a cohort of individuals (see, for example, [11,12,13,14,15,16,17]). We have also developed a number of physics-informed (based on the first principles diffusion-based model given in (1.1)–(1.5) below) data-based machine learning schemes using hidden Markov models, generative adversarial and long short-term memory neural network models for estimating BrAC from biosensor TAC data, which are also trained using population data and do not require individual calibration [18,19,20].

    In all of our approaches to the TAC to BAC/BrAC conversion problem, the underlying model was taken to be based on the first principles physics based initial-boundary value problem for a parabolic partial differential equation. This will also be the basic model to which we will direct our efforts in the present treatment. Transformed to be in terms of dimensionless variables, the model is given by

    xt(t,η)=q12xη2(t,η),0<η<1,t>0, (1.1)
    q1xη(t,0)=x(t,0),t>0, (1.2)
    q1xη(t,1)=q2u(t),t>0, (1.3)
    x(0,η)=x0,0<η<1, (1.4)
    y(t)=x(t,0),t>0, (1.5)

    where t and η are the temporal and spatial variables, respectively, and x(t,η) indicates the concentration of ethanol in the epidermal layer of the skin at time t and depth η, where η=0 is at the skin surface and η=1 is at the boundary between the epidermal and dermal layers of the skin. The input to the system is u(t), which is the BrAC/BAC at time t, and the output is y(t), which is the TAC at time t. Equation (1.1) represents the transport of ethanol through the epidermal layer of the skin. The boundary conditions (1.2) and (1.3) represent respectively the evaporation of ethanol at the skin surface and the flux of ethanol across the boundary between the epidermal and dermal layers. It is assumed that there is no alcohol in the epidermal layer of the skin at time t=0, so the initial condition (1.4) is x0(η)=0, 0<η<1. Finally, the output Eq (1.5) represents the TAC level measured by the biosensor at the skin surface.

    The parameters in the system (1.1)–(1.5) that will be assumed to be random are q1 and q2, which represent respectively the normalized diffusivity and the normalized flux gain at the boundary between the dermal and epidermal layers. The values or distributions of these parameters are assumed to depend on environmental conditions, the particular sensor being used, and the physiological characteristics of the individual wearing the sensor. The parameter vector is q=(q1,q2)Q, where Q is assumed to be a compact subset of R+×R+ with metric d.

    In population modeling, we can statistically classify the methods as parametric or nonparametric. In the parametric approach, we assume that the general structure of the distribution is known a-priori but with unknown parameters. Then, for example, if we know that the distribution is normal with unknown mean and variance, the estimation problem is to estimate these two unknown parameters. On the other hand; in the nonparametric approach, the structure of the distribution is assumed to be unknown, and the problem is to estimate the distribution itself. In either of these paradigms, different statistical approaches to the estimation problem can be taken. For example, in [15,16,17], a parametric least squares naive pooled data approach was used, while in [21,22], the approach was nonparametric. In [12], a Bayesian framework was developed, and in the present treatment we consider a mixed effects (see, for example, [23,24,25]) maximum likelihood-based statistical model. In the mixed effects model, it is assumed that observations are specific to a single individual plus a random error. The mixed effects model is a combination of the fixed-effects model, which describes the characteristics for an average individual in the population, and the random-effects model, which describes the inter-individual variability [26]. An overview of these different statistical approaches in the context of pharmacokinetics can be found in [27].

    While TAC and BrAC are related, they are not precisely equivalent in terms of quantification. BrAC serves as a reliable indicator of BAC, offering valuable quantitative insights into the immediate impact of alcohol on judgment, motor skills and cognition. On the contrary, TAC lags behind BrAC because it must pass through the skin before measurement. TAC levels exhibit variability among individuals, devices and across drinking episodes, in contrast to the relative consistency of BrAC across different people and conditions. Consequently, interpreting TAC in relation to real-time blood alcohol levels becomes challenging, limiting its usefulness as a quantitative alcohol measure. Therefore, having a method to estimate BrAC from TAC would enhance its practicality for researchers, healthcare professionals and individuals seeking meaningful and time-sensitive information about their alcohol levels. The inverse problem of estimating the input BrAC signal u to the system (1.1)–(1.5) based on observations of the TAC output signal y is the central focus of our team's research effort. In our treatment here, we consider only the problem of estimating the distribution of the random parameters in the population model and how the uncertainty in the parameters propagates through the system of Eqs (1.1)–(1.5) to the TAC output y. Once the distributions of the random parameters in the the model (1.1)–(1.5) have been estimated, we have previously developed (and continue to investigate) a number of deconvolution or inverse filtering methods that can then be used to obtain an estimate for the BrAC signal along with a credible or error band that result from the uncertainty in the model parameters as described by their estimated distributions. We have looked at a number of approaches including frequency domain, Bayesian, linear-quadratic Gaussian compensator, maximum likelihood, autoregressive moving average (ARMA) and nonlinear least squares based techniques [11,28,29,30].

    In addition to the work of our group on TAC to BAC/BrAC conversion cited above, other researchers have also been looking at this problem and have tried a number of different approaches. For example, in [31,32], a more traditional approach based on standard linear regression techniques was developed and discussed. A number of ideas from the machine learning literature have also been considered. In [33], a scheme based on random forests was used to recover BrAC from TAC, and in our group, there are the machine learning based treatments cited earlier [18,19,20].

    An outline of the remainder of the paper is as follows. In Section 2, we provide a summary of the Prohorov metric on the set of probability measures as it was used by Banks and his coauthors in [21]. In Section 3, we define our mathematical model in the form of a random discrete-time dynamical system and we define the maximum likelihood estimator for the distribution of the random parameters. In Section 4, we establish the existence and consistency of the maximum likelihood estimator for the distribution of the random parameters, while in Section 5, we demonstrate the convergence of finite dimensional approximations for our estimator. In Section 6, we summarize results for abstract parabolic systems, their finite dimensional approximation and an associated convergence theory. In Section 7, the application of our scheme to the transdermal transport of alcohol is presented and discussed. This includes numerical studies for two examples, one involving simulated data and the other actual data collected in the laboratory of one of the coauthors, Dr. Susan Luczak in the Department of Psychology at University of Southern California (USC). For the simulated data example in Section 7.1, we are able to observe the convergence of the estimated distribution of the random parameter vector q=(q1,q2) to the "true" distribution as the number of drinking episodes increases, as the number of Dirac measure nodes increases and as the level of discretization in the finite dimensional approximations increases. In the actual data example discussed in Section 7.2, we apply the leave-one-out cross-validation (LOOCV) method by first estimating the distribution of the parameter vector q using a training set, and then estimating the TAC output using the estimated distribution and the BrAC input of a testing episode.

    Banks and his coauthors developed a framework for estimation of the probability measure for random parameters in continuous-time dynamical systems based on the Prohorov metric [21]. Here, we summarize the Prohorov metric and its properties.

    Let Q be a Hausdorff metric space with metric d. Define

    Cb(Q)={f:QR|fis bounded and continuous},

    and given any probability measure PP(Q), where P(Q) denotes the set of all probability measures defined on ΣQ, the Borel sigma algebra on Q, and some ε>0, an ε-neighborhood of P is defined by

    Bε(P)={˜P||Qf(q)d˜P(q)Qf(q)dP(q)|<ε,for allfCb(Q)}.

    Let EΣQ, and define the ε-neighborhood of E by

    Eε={˜qQ|d(˜q,E)<ε}={˜qQ|infqEd(q,˜q)<ε}.

    Given two probability measures, P and ˜P in P(Q), the Prohorov metric ρ on P(Q)×P(Q) is defined such that

    ˜PBε(P)ρ(P,˜P)<ε,

    where

    ρ(P,˜P)=inf{ε>0|˜P(E)P(Eε)+εandP(E)˜P(Eε)+ε,for allEΣQ}.

    It can be shown that (P(Q),ρ) is a metric space. Also, the Prohorov metric metrizes the weak convergence of measures, i.e., given a sequence of measures PMP(Q), for all M=1,2,, and PP(Q),

    PMwPρ(PM,P)0.

    It is important to note that the weak topology and the weak topology are equivalent on the space of probability measures.

    For some nq, QRnq and PP(Q), consider the random vector X:QRnq on the probability space (Q,ΣQ,P) given by X(q)=q for qQ. The cumulative distribution function for X is given by FX(q1,,qn)=P(X×nq=1(,q])=P(×nq=1{(,q]Q}). In this case, it follows that if PM,P0(P(Q), ρ), for M=1,2,, then ρ(PM,P0)0 if and only if FXMFX0 at all points of continuity of FX0. Consequently, Prohorov metric convergence and weak and weak convergence in P(Q) are also referred to as convergence in distribution.

    If q1,q2Q, then ρ(δq1,δq2)=min{d(q1,q2),1}, where δqjD={δq|qQ}, the space of Dirac measures on Q, where for all EΣQ,

    δq(E)={1if qE0if qE

    The metric space (P(Q),ρ) is separable if and only if the metric space (Q,d) is separable. The sequence {qj}j=1 is Cauchy in (Q,d) if and only if the sequence {δqj}j=1 is Cauchy in (P(Q),ρ). We also have that (Q,d) is complete if and only if (P(Q),ρ) is complete, and (Q,d) is compact if and only if (P(Q),ρ) is compact. The details and proofs can be found in [21].

    Assume the metric space (Q,d) is separable and let Qd={qj}j=1 be a countable dense subset of Q. Define the dense (see [21]) subset of P(Q), ˜Pd(Q), as

    ˜Pd(Q)={PP(Q)|P=Mj=1pjδqj,qjQd,MN,pj[0,1]Q,Mj=1pj=1}, (2.1)

    the collection of all convex combinations of Dirac measures on Q with rational weights pj at nodes qjQd, and for each MN let

    PM(Q)={P˜Pd(Q)|P=Mj=1pjδqj,qj{qj}Mj=1}.

    Consider the following discrete-time mathematical model for the ith subject at time-step k

    xk,i(qi)=gk1(xk1,i(qi),uk1,i;qi),k=1,,ni,i=1,,m,x0,i=ϕ0,i,i=1,,m,

    where qi is the ith subject's parameter vector in Q, denoting the set of admissible parameters, gk1:H×Rν×QH, H is in general an infinite dimensional Hilbert space, and uk1,iRν is the input. The output is given by

    yk,i(qi)=hk(xk,i(qi),ϕ0,i,uk,i;qi),k=1,,ni,i=1,,m,

    where hk:H×H×Rν×QR.

    For the mixed effects model, we define

    Yk,i=yk,i(qi)+ek,i,k=1,,ni,i=1,,m, (3.1)

    where for each i=1,2,,m, ek,i are independent and identically distributed (i.i.d.) with mean zero, variance σ2 and ek,iφ, k=1,2,,ni, where φ is a density with respect to a sigma finite measure μ on R and assumed to be continuous on R. We assume that the random vectors [e1,i,,eni,i] are independent with respect to i, i=1,2,,m; that is, the error is independent across individuals and conditionally independent within individuals (i.e., given qi). For each i=1,2,,m, let Yi=[Y1,i,Y2,i,,Yni,i]T, yi(qi)=[y1,i(qi),y2,i(qi),,yni,i(qi)]T, ei=[e1,i,e2,i,,eni,i]T and rewrite (3.1) as

    Yi=yi(qi)+ei,i=1,,m.

    Then, for i=1,2,,m, Yi are independent with Yifi(;qi):RniR, where

    fi(v;qi)=nik=1φ(vkyk,i(qi)),i=1,,m, (3.2)

    where v=[v1,v2,,vni]TRni. Let PP(Q) denote a probability measure on ΣQ, where P(Q) denotes the set of all probability measures defined on ΣQ, and let P0P(Q) be the "true" distribution of the random vector qi. The goal is to find an estimate of P0. In order to generate an estimator for P0 and establish theoretical results and computational tools, we use the nonparametric maximum likelihood (NPML) approach introduced by Lindsay and Mallet in [34,35], as well as the Prohorov metric-based framework on P(Q) introduced by Banks and his coauthors in [21], summarized in Section 2.

    For PP(Q) and i=1,2,,m, let

    Li(P;Yi)=Qfi(Yi;qi)dP(qi)=Qnik=1φ(Yk,ihk(xk,i(qi),ϕ0,i,uk,i;qi))dP(qi) (3.3)

    be the contribution of the ith subject to the likelihood function

    Ln,m(P;Y)=mi=1Li(P;Yi)=mi=1Qfi(Yi;qi)dP(qi)=mi=1Qnik=1φ(Yk,ihk(xk,i(qi),ϕ0,i,uk,i;qi))dP(qi), (3.4)

    where n={ni}mi=1 and Y={Yi}mi=1. The goal is to find P that maximizes the likelihood function.

    Define the estimator

    Pn,m=argmaxPP(Q)Ln,m(P;Y)=argmaxPP(Q)mi=1Qfi(Yi;qi)dP(qi)=argmaxPP(Q)mi=1Qnik=1φ(Yk,ihk(xk,i(qi),ϕ0,i,uk,i;qi))dP(qi). (3.5)

    Let ^𝓎k,i be realizations of the random variables Yk,i, and define

    (3.6)

    where , with .

    The results of Lindsay and Mallet in [34,35] states that the maximum likelihood estimator ˆPn,m can be found in the class of discrete distributions with at most m support points; that is, ˆPn,mPM(Q), where Mm. So, we define our approximating estimator over the set PM(Q) where M denotes the number of nodes, {qj}Mj=1. As a result, the optimization is over a finite set of parameters, being the rational weights {pj}Mj=1.

    Also, we cannot exactly compute the maximum likelihood estimator ˆPn,m since yk,i must be approximated numerically by yNk,i using a Galerkin numerical scheme with N denoting the level of discretization. Thus, our approximating estimator is

    (3.7)

    In [34], the existence and uniqueness of a maximum likelihood estimator of a mixing distribution using the geometry of mixture likelihoods was established. Similarly, in [35], the existence and uniqueness of the maximum likelihood estimator for the distribution of the parameters of a random coefficient regression model was established. Here, we provide an existence argument based on the maximization of a continuous function over a compact set.

    The following theorem establishes the existence of the estimator ˆPn,m in (3.7), obtained from the realizations {^𝓎k,i},k=1,,ni,i=1,,m of the random variables {Yk,i},k=1,,ni,i=1,,m. This is sufficient for establishing the existence of the maximum likelihood estimator Pn,m in (3.5).

    Theorem 4.1. For i=1,2,,m, let Li be given by Eq (3.3) and let be given by equation (3.4), where with . Assume that for each , we have a continuous function , and also for each PP(Q) we have a measurable function Li(P;.):RniR. Then there exists a measurable function ˆPn,m:mi=1RniP(Q) such that

    Proof. The theorem can be proven in a similar way as in [21] with the difference that we are taking the sup (instead of the inf) of a continuous function over a compact set.

    In order to establish the consistency of the maximum likelihood estimator Pn,m, we show that ρ(Pn,m,P0) converges almost surely to zero. We do this by applying a theorem by Kiefer and Wolfowitz in [36], establishing that the nonparametric maximum likelihood approach is statistically consistent. In other words, as the number of subjects m gets larger, the estimator Pn,m converges in probability to P0, the "true" distribution in the sense of the Prohorov metric, or weakly or in distribution. Here, we have set up our problem in a way that makes establishing the consistency a straightforward application of the consistency result in [36].

    Theorem 4.2. For each i=1,...,m, assume that the map qifi(v;qi) from Q into R is continuous for each vRni, and fi(v;qi) is measurable in v for any qiQ, where fi is given by (3.2). Assume further that P0 is identifiable; that is, for P1P(Q) with P1P0, we have

    mi=1zi0Qnik=1φ(zkhk(xk,i(qi),ϕ0,i,uk,i;qi))dP1(qi)dμnimi=1zi0Qnik=1φ(zkhk(xk,i(qi),ϕ0,i,uk,i;qi))dP0(qi)dμni,

    for at least one z=[zT1,,zTm]TRmi=1ni, where for i=1,2,,m, zi=[z1,,zni]TRni, and the technical integrability assumption holds; that is, for any PP(Q),

    limε0EP0[logsup˜PBε(P)Li(˜P;Yi)Li(P0;Yi)]+<,

    where Li is given by Eq (3.3). Then, as m, ρ(Pn,m,P0)0 almost surely (i.e., with probability 1), and therefore in probability as well.

    Proof. The assumptions we have made in the previous section and in the statement of the theorem are sufficient to argue that Assumptions 1–5 in [36] are satisfied. The conclusion of the consistency result in [36] is that the cumulative distribution functions Fn,m corresponding to Pn,m converge almost surely to the cumulative distribution function F0 corresponding to P0 at every point of continuity of F0. It follows that ρ(Pn,m,P0)0 almost surely (i.e., with probability 1); thus, in probability as well, as m, and the theorem is proven.

    We want to establish the convergence of the finite dimensional maximum likelihood estimators to the maximum likelihood estimator corresponding to the infinite dimensional model. As mentioned earlier, we cannot actually compute ˆPn,m in (3.6) and, consequently, we approximate it by ˆPNn,m,M in (3.7). Consider the following assumptions:

    A1. For all n, m, and N, the map is a continuous map.

    A2. For any PM,PP(Q), M=1,2,, such that ρ(PM,P)0, we have as N,M for i=1,,m.

    A3. For all PP(Q) and i=1,,m, and are uniformly bounded.

    Theorem 5.1. Under assumptions A1–A3, there exists maximizers ˆPNn,m,M given by Eq (3.7). In addition, there exists a subsequence of ˆPNn,m,M that converges to ˆPn,m given by Eq (3.6) as M,N.

    Proof. For all n, m, and N, by continuity of the map per assumption A1 and compactness of (P(Q),ρ), we can conclude that ˆPNn,m,M exists.

    In [21], it is shown that ˜Pd(Q) given by Eq (2.1) is a dense subset of P(Q). Thus, for M=1,2,, construct a sequence of probability measures PMPM(Q)˜Pd(Q)P(Q), such that ρ(PM,P)0 in P(Q). Then, by assumptions A2 and A3, we have

    Consequently, as M,N.

    In addition, by definition, for each n, m and N and for all PMPM(Q), we have

    (5.1)

    In addition, by compactness of P(Q), there exists a subsequence of ˆPNn,m,M that converges to ˆPn,m as M,N. Thus, by taking the limit in (5.1) as M,N, for all PP(Q), we find that

    thus, as given in Eq (3.6).

    In practice, to achieve a desired level of accuracy, M and N are fixed sufficiently large. We choose a sufficiently large value for N. How large that needs to be, of course, depends on the particular numerical discretization scheme chosen. The most common choice would be using a Galerkin-based method to the approximate yk,i by yNk,i, where

    yNk,i=hk(xNk,i(qi),ϕN0,i,uk,i;qi),

    which denotes the discretization of the output of the model.

    We also choose a sufficiently large value for M, the number of nodes, {qj}Mj=1. Therefore, the optimization problem is reduced to a standard constrained estimation problem over Euclidean M-space, in which we determine the values of the weights pj at each node qj with the constraints that they all be nonnegative and sum to one. By Eq (3.7) it follows that

    where ˜p=(p1,,pM)~RM={˜p|pjR+,Mj=1pj=1}.

    We note that computing ˆPNn,m,M involves high order products of very small numbers that, not unexpectedly, can cause numerical underflow. In order to mitigate this, we maximize the log-likelihood function instead and rewrite it in a form that lends itself to the use of the MATLAB optimization routine logsumexp as follows:

    ˆPNn,m,M=argmax˜p~RMlog(mi=1Mj=1pjnik=1φ(^𝓎k,ihk(xNk,i,ϕN0,i,uk,i;qj)))=argmax˜p~RMmi=1log(Mj=1pjnik=1φ(^𝓎k,ihk(xNk,i,ϕN0,i,uk,i;qj)))=argmax˜p~RMmi=1log(Mj=1exp(log(pjnik=1φ(^𝓎k,ihk(xNk,i,ϕN0,i,uk,i;qj)))))=argmax˜p~RMmi=1log(Mj=1exp(log(pj)+nik=1log(φ(^𝓎k,ihk(xNk,i,ϕN0,i,uk,i;qj))))). (5.2)

    In order to apply our estimation theory to Eqs (1.1)–(1.5), our model for the transdermal transport of ethanol given in Section 1, we reformulate it as an abstract parabolic system. We briefly describe what an abstract parabolic system is, its properties and its finite dimensional approximation, and then we show how assumptions A1A3 are satisfied for such a system.

    Let H and V be Hilbert spaces with V densely and continuously embedded in H. Pivoting on H, it follows that H is therefore densely and continuously embedded in the dual of V, V. This is known as a Gelfand triple and is generally written as VHV [37]. An abstract parabolic system is a dynamical system of the following form

    <˙x,ψ>V,V+a(q;x,ψ)=<B(q)u,ψ>V,V,ψV,x(0)=x0,y(t)=C(q)x(t), (6.1)

    where <,>V,V denotes the duality pairing between V and V, Q is as defined in Section 1, and for each qQ, a(q;.,.):V×VC is a sesquilinear form satisfying the following three assumptions:

    B1. (Boundedness) There exists a constant α0 such that for all ψ1,ψ2V, we have

    |a(q;ψ1,ψ2)|α0ψ1Vψ2V.

    B2. (Coercivity) There exists λ0R and μ0>0 such that for all ψV, we have

    a(q;ψ,ψ)+λ0|ψ|2Hμ0ψ2V.

    B3. (Continuity) For all ψ1,ψ2V and q,˜qQ, we have

    |a(q,ψ1,ψ2)a(˜q,ψ1,ψ2)|d(q,˜q)ψ1Vψ2V.

    In these assumptions, .V and |.|H denote the norm on the spaces V and H, respectively. Further, in (6.1), B(q):RνV, and C(q):VR are bounded linear operators with initial conditions x0H, input uL2([0,T],Rν), and output yL2([0,T],R).

    It can be shown that the system in (6.1) has a unique solution in

    {ψ|ψL2([0,T],V),˙ψL2([0,T],V)}C([0,T],H)

    using standard variational arguments (such as in [38]). However, we use a linear semigroup approach to convert the system in (6.1) into a discrete-time state space model, and then use arguments from linear semigroup theory [6,39] to argue convergence of finite dimensional Galerkin-based approximations and conclude that assumptions A1A3 are satisfied.

    Assumptions B1 and B2 yield that the form a(q;.,.) defines a bounded linear operator A(q):VV given by

    <A(q)ψ1,ψ2>V,V=a(q;ψ1,ψ2),

    for ψ1,ψ2V. If we restrict the operator A(q) to the subspace

    Dom(A(q))={ψV|A(q)ψH},

    it becomes the infinitesimal generator of a holomorphic or analytic semigroup {eA(q)t|t0} of bounded linear operators on H. The operator A(q) is referred to as being regularly dissipative [5,6,37]. Moreover, this semigroup can also be extended and restricted to be a holomorphic semigroup on V and V, respectively [5,37].

    The system in (6.1) can now be written in state space form with time invariant operators A(q), B(q) and C(q) as

    ˙x(t)=A(q)x(t)+B(q)u(t),x(0)=x0,y(t)=C(q)x(t). (6.2)

    The operator form of the variation of constants formula then yields what is known as a mild solution of (6.2) given by

    x(t;q)=eA(q)tx0+t0eA(q)(ts)B(q)u(s)ds,t0,y(t;q)=C(q)x(t;q). (6.3)

    To obtain the corresponding discrete or sampled time form of the system given in (6.2) or (6.3), let τ>0 be the length of the sampling interval and consider strictly zero-order hold inputs of the form u(t)=uk1,t[(k1)τ,kτ),k=1,2,. Then, let xk=x(kτ) and yk=y(kτ),k=1,2,. By applying (6.3) on each subinterval [(k1)τ,kτ], k=1,2,, we obtain the discrete-time dynamical system given by

    xk=ˆA(q)xk1+ˆB(q)uk1,k=1,2,, (6.4)
    yk=ˆC(q)xk,k=1,2,, (6.5)

    where x0V, ˆA(q)=eA(q)τ, ˆB(q)=τ0eA(q)sB(q)ds and ˆC(q)=C(q).

    Using a standard Galerkin approach [40], we can approximate the discrete-time system given in (6.4)–(6.5) by a sequence of approximating finite dimensional discrete-time systems in a sequence of finite dimensional subspaces VN of V. In order to argue convergence, we will require the following additional assumption concerning the subspaces VN.

    C1. (Approximation) For every xV, there exists xNVN such that xxNV0 as N.

    We consider the sequence of approximating finite dimensional discrete-time systems by

    xNk=ˆAN(q)xNk1+ˆBN(q)uk1,k=1,2,,yNk=ˆCN(q)xNk,k=1,2,,

    where ˆAN(q)=eAN(q)τ, ˆBN(q)=τ0eAN(q)sBN(q)ds and ˆCN(q)xk=ˆC(q), where for each qQ, AN(q) is the linear operator on VN obtained by restricting the form a(q;.,.) to VN×VN, i.e., for ψN1,ψN2VN,

    <AN(q)ψN1,ψN2>V,V=a(q;ψN1,ψN2).

    In addition, BN(q)=πNB(q), where in this definition, πN is the natural extension of the orthogonal projection operator πN:HVN to V from its dense subspace H. We also set xN0=πNx0VN.

    Under the assumptions B1–B3 and C1 using the Trotter-Kato approximation theorem from the theory of linear semigroups of operators [39,41], we were able to conclude that limNxNkxkV=0 and limN|yNkyk|=0 for each x0V, and uniformly in q for qQ and k{1,2,,K}, for any fixed KN+.

    We can now use the results described in the previous paragraphs to show that an abstract parabolic system satisfies assumptions A1–A3 given in Section 5. To show that the assumption A1 is satisfied, we need to show that for all n, m, and N, the map is a continuous map. It suffices to show that for any fixed n, m, and N, and for any sequence of probability measures PM, such that ρ(PM,P)0 in P(Q), we have as M. Toward this end, we see that

    by definition of the Prohorov metric. It follows that the assumption A1 is satisfied.

    Next, we show that the assumption A2 is satisfied. We have that limNxNk,ixk,iV=0 and limN|yNk,iyk,i|=0 for each x0,iV, uniformly in qi for qiQ, k=1,,ni,i=1,,m. We want to show that for any sequence of probability measures PM, such that ρ(PM,P)0 in P(Q), and for i=1,,m, as M,N, we have .

    Recall that φ is assumed to be continuous. Let ε>0. Choose N0 such that for NN0, and for every M, |Q(nik=1φ(^𝓎k,iyNk,i(qi))nik=1φ(^𝓎k,iyk,i(qi)))dPM(qi)|<ε/2. Then, we have

    where the second term is less than ε/2 by definition of the Prohorov metric. Consequently, the assumption A2 is satisfied.

    Finally, we want to show that the assumption A3 is satisfied. We want to show that for all PP(Q) and for i=1,,m, and are uniformly bounded. Recall that the parameter space Q is compact. Thus, for qiQ, and for each N, yNk,i(qi) are uniformly bounded. Similarly, yk,i(qi) are also uniformly bounded and we also have that |yNk,i(qi)yk,i(qi)|0 uniformly in qi for qiQ. Therefore, we can conclude that the assumption A3 is satisfied.

    To apply the results established in Section 6 to the system (1.1)–(1.5) in Section 1, the system must first be written in weak form, then the parameter space Q, the Hilbert spaces H and V, the sesquilinear form a(q;.,.) and the operators B(q) and C(q) must all be identified. Also, the approximating subspaces VN must be chosen, and finally assumptions B1–B3 and C1 must all be shown to be satisfied.

    The parameter space Q is assumed to be a compact subset of R+×R+ with any p-metric denoted by dQ. Let x=(x1,,xn)Rn and y=(y1,,yn)Rn, then the p-metric is defined by

    dp(x,y)=(ni=1|xiyi|p)1/p

    for any p[1,). For p=, the p-metric is defined by

    dp(x,y)=maxi=1,,n|xiyi|.

    Let H=L2(0,1) and V=H1(0,1) with their standard inner products and norms. It follows that V=H1(0,1) and the three spaces H, V and V form a Gelfand triple. To rewrite the system (1.1)–(1.5) in weak form, we multiply by a test function ψV and integrate by parts to obtain

    <˙x(t),ψ>V,V+10q1xη(t,η)ψ(η)dη+x(t,0)ψ(t,0)=q2u(t)ψ(1),

    where <,>V,V denotes the duality pairing between V and V. Then, for qQ, uR and ˜ψ,ψV, we set

    a(q;˜ψ,ψ)=10q1˜ψ(η)ψ(η)dη+˜ψ(0)ψ(0),<B(q)u,ψ>V,V=q2uψ(1),C(q)ψ=Cψ=ψ(0).

    We can establish that assumptions B1–B3 are satisfied using arguments involving the Sobolev embedding theorem (see [42]). Also, the operators B(q) and C(q) are continuous in the uniform operator topology with respect to qQ. It follows from Section 6 that

    gk1(xk1,i(qi),uk1,i;qi)=ˆA(qi)xk1,i(qi)+ˆB(qi)uk1,i,k=1,2,,ni,i=1,,m,hk(xk,i(qi),ϕ0,i,uk,i;qi)=ˆC(qi)xk,i(qi),k=1,,ni,i=1,,m,

    where ˆA(q)=eA(q)τ, ˆB(q)=τ0eA(q)sB(q)ds and ˆC(q)=C(q), with τ>0 which is the length of the sampling interval.

    Let VN,N=1,2,, be the span of the standard linear splines defined with respect to the uniform mesh {0,1/N,2/N,,(N1)/N,1} on [0,1]. Then, assumption C1 is satisfied by standard arguments for spline functions (see, for example, [43]). If for each i=1,2,,m, we define xNk,i and yNk,i as in (6.4)–(6.5), then by the arguments at the end of Section 6, we conclude that assumptions A1–A3 are satisfied.

    In the following two subsections, 7.1 and 7.2, we present the application of our scheme to the transdermal transport of alcohol in two examples, one involving simulated data, and the other using actual human subject data collected in the Luczak laboratory at USC. For the simulated data, we want to show the convergence of the estimated distribution of the parameter vector q=(q1,q2) to the "true" distribution as the number of drinking episodes increases, as the number of nodes increases and as the level of discretization in the finite dimensional approximations increases. For the actual data we apply the LOOCV method by first estimating the distribution of the parameter vector q using the training set, and then we estimated the TAC output using the estimated distribution of the parameter vector q and the BrAC input of the test set.

    In this example, we estimate the distribution of the parameter vector q=(q1,q2) in the system (1.1)–(1.5) by first simulating TAC data in MATLAB with the assumption that the two parameters q1 and q2 are i.i.d. with a Beta distribution, q1,q2Beta(2,5). Thus, their joint cumulative distribution function (cdf) is the product of their marginal Beta(2,5) cdfs.

    From Eq (3.1), we have

    Yk,i=yk,i(qi)+ek,i,k=1,,ni,i=1,,m,

    where m is the number of drinking episodes and yk,i(qi) is the observed TAC for the ith drinking episode at time step k. We let P0 be the product of the cdfs of two independent Beta(2,5) distributions, and ek,iN(0,106),k=1,,ni,i=1,,m. The reason for the small variance is because the noise accounts for errors in the measurements caused by factors that result in small perturbations to the measurements, such as humidity or temperature as external environmental factors or the level of skin moisture or greasiness among different drinking episodes. It is assumed that these types of factors are normally distributed, causing small variations to the measurements.

    To approximate the PDE model for the TAC observations, we used the linear spline-based Galerkin approximation scheme described in Section 6 with N equally spaced subintervals from [0,1] (see [15,16,17]). We want to compute ˆPNn,m,M given by Eq (5.2), where qj=(qj1,qj2) is chosen as M uniform meshgrid coordinates on [0,1]×[0,1]. We make the assumption that there is no alcohol in the epidermal layer of the skin at time t=0, so we let ϕN0,i=0. The constrained optimization problem over Euclidean M-space was solved using constrained optimization routine FMINCON (find minimum of constrained nonlinear multivariable function) from the Optimization Toolbox in MATLAB applied to the negative of the log-likelihood function.

    We note that before turning to a mixed effects statistical model, we had considered an output or observation that was aggregated TAC; in this case, the appropriate underlying statistical model was the naive pooled model (i.e., the data point for each drinking episode at a certain time is an observation of the mean behavior plus a random error). However, it is not difficult to show that when this observation is used, if the q1 and q2 are assumed to be independent, then their joint distribution is not identifiable. Consequently, when the nonlinear least squares-based constrained optimization problem is solved, the inherent ill-posedness of the inverse problem results in undesirable oscillations. To mitigate this behavior, we had to introduce an appropriately weighted regularization term in the performance index being minimized. One advantage of the mixed effects statistical model presented here is that regularization is not required.

    In order to demonstrate the consistency of our estimator, we show that as m, the number of drinking episodes, increases, the estimated cdf of the parameter vector q=(q1,q2) approaches the "true" cdf, the product of two Beta(2,5) cdfs. In order to simulate realistic longitudinal TAC vectors representing data that might be collected by the TAC biosensor for an individual's drinking episode, we used BrAC data collected in the Luczak laboratory as the input to the model and generated random samples of q1 and q2, i.i.d. Beta(2,5) in MATLAB. Using the algorithm developed in the current paper, we estimated the distribution of the random parameter vector by solving the optimization problem for different cases based on the number of drinking episodes m and observed the convergence of the estimated distribution to the "true" distribution as m increases.

    To quantify this, let D be the sum of the squared differences at each node between the estimated and the "true" distribution, the product of two Beta(2,5) cdfs. Let pj and bj be the weights at the node qj of the estimated and "true" distribution, respectively. So, we have

    D=Mj=1(pjbj)2

    representing the error in estimating the weights at each node. This is a metric being used to quantify that the difference between the estimated distribution and the "true" distribution decreases as the number of drinking episodes increases.

    We fixed the number of nodes M and the level of discretization N sufficiently large. We set M=400 and N=128. We estimated the distribution of the parameter vector for different cases based on different numbers of drinking episodes, m{1,3,7,9,16,42}, and calculated D for each case. Our results are summarized in Table 1. We observed that as the number of drinking episodes m increases, the sum of the squared differences at each node between the estimated distribution and the "true" distribution D decreases.

    Table 1.  Decrease in D, the sum of the squared differences at each node between the estimated and the "true" distribution, with increasing m, the number of drinking episodes, for fixed values of the number of nodes M=400 and the level of discretization N=128.
    m D
    1 39.0164
    3 28.3091
    7 8.3247
    9 7.1750
    16 3.5697
    42 3.0337

     | Show Table
    DownLoad: CSV

    In Figure 2, each row of the figure contains three different views of the same plot of the estimated distribution and the "true" distribution (again, the product of two Beta(2,5) cdfs) for different numbers of drinking episodes m=7, m=16 and m=42 in the top, middle and bottom rows, respectively, with the number of nodes set to M=400 and the level of discretization to N=128. We observe that as m increases, our estimated distribution gets "closer" to the "true" distribution, which agrees with the numerical results that are shown in Table 1.

    Figure 2.  Each row of the figure contains three different views of the same plot of the estimated distribution and the "true" joint Beta(2,5) distribution for different numbers of drinking episodes m=7, m=16 and m=42 in the top, middle and bottom rows, respectively, for a fixed number of nodes M=400 and level of discretization N=128. We observe that as m increases from the top row to the bottom row, our estimated distribution gets "closer" to the "true" distribution, which agrees with the numerical results that are shown in Table 1.

    Next, we show that as the number of nodes M and the level of discretization N increases, the normalized sum of squared differences at each node between the estimated and the "true" distribution decreases.

    First, we fixed the level of discretization at N=128 and we increased the number of nodes M. Let

    ˉDM=1MMj=1(pjbj)2.

    In Table 2, we observe that for the fixed value of N=128, as the number of nodes M increases, the normalized sum of the squared differences at each node between the estimated distribution and the "true" distribution ˉDM decreases.

    Table 2.  Decrease in ˉDM, the normalized sum of the squared differences at each node between the estimated distribution and the "true" distribution, with increasing M, the number of nodes, for a fixed value for the level of discretization N=128.
    M ˉDM
    25 0.03795
    100 0.03025
    225 0.02974
    400 0.02081

     | Show Table
    DownLoad: CSV

    Next, we fixed the number of nodes at M=400 and we increased the level of discretization N. Let

    ˉDN=1NMj=1(pjbj)2.

    In Table 3, we observe that for N=128 fixed, as the number of nodes M increases, the normalized sum of the squared differences at each node between the estimated distribution and the "true" distribution ˉDN decreases.

    Table 3.  Decrease in ˉDN, the normalized sum of the squared differences at each node between the estimated distribution and the "true" distribution, with increasing N, the level of discretization, with the number of nodes fixed at M=400.
    N ˉDN
    4 1.22783
    16 0.65644
    64 0.18565
    128 0.06504

     | Show Table
    DownLoad: CSV

    The choice of the "true" distribution for q1 and q2 in the simulation case is the scatterplot of samples for a set of 18 drinking episodes, including BrAC and TAC measurements of different individuals obtained by a deterministic approach in [22]. However, the Beta(2,5) distribution was chosen strictly for the purpose of demonstration. When applying our algorithm to actual clinic or lab collected human subject data, a valuable feature of our developed methodology in using a nonparametric approach is that we do not need to make any assumptions about the family of feasible distributions for the parameter vector unlike the parametric approach in [15,16,17]. In addition, the independent and identically distributed assumption was also very simplistic given that q1 and q2 parameters depend on the same individual and environmental conditions at the time of measurements. In the nonparametric approach, there is no need to make any assumptions about the distribution, resulting in relaxing the restrictive assumption of the distribution being independent and identically distributed. Thus, this assumption is relaxed in the flexible nonparametric approach used in the next example.

    The two datasets used in this example were obtained by two different alcohol biosensors: SCRAM CAM and WrisTASTM7 [44,45]. We fixed the number of nodes at M=400 and the level of discretization at N=128, both sufficiently large with respect to convergence as we observed in our simulation data examples. From each dataset, we chose m=9 different drinking episodes. We split the drinking episodes into a training set consisting of eight drinking episodes and a testing set consisting of one drinking episode. This way, we could apply the LOOCV method. We repeated this partitioning process nine times, each time leaving out a different drinking episode. Using the training set, we first estimated the distribution of the parameter vector q=(q1,q2). Next, we sampled 100 parameter vectors q=(q1,q2) from the estimated distribution, and using those along with the BrAC input from the testing dataset, we simulated 100 TAC longitudinal signals. From these 100 simulated TAC signals, we estimated the "true" TAC by computing the mean at each time, and we provided what we refer to as a 95 conservative error band, or simply as a 95 error band, by taking the 2.5 and 97.5 percentiles. This approach for the error band is also used for a number of statistics associated with the TAC curve that are of particular interest to researchers and clinicians working in the area of alcohol use disorder.

    For the first example, we considered the dataset collected using the SCRAM alcohol biosensor. Prior to applying the LOOCV method, in order to visualize the estimated density and distribution of q=(q1,q2) and the marginal densities of q1 and q2, we trained the algorithm on all nine drinking episodes. Figure 3 illustrates the four aforementioned plots. In the estimated density plot, we can see that our numerical result for this example is in agreement with the theoretical result in Lindsay and Mallet [34,35], which states that the maximum likelihood estimator ˆPn,m can be found in the class of discrete distributions with at most m support points, i.e., ˆPn,mPM(Q), where Mm. In this example, since we had nine drinking episodes, the estimated density plot displays the support points among 400 nodes for q=(q1,q2).

    Figure 3.  The estimated probability density function (pdf) (top left), the estimated cdf (top right), marginal density of q1 (bottom left), and marginal density of q2 (bottom right) obtained from m=9 drinking episodes collected using the SCRAM alcohol biosensor, for a fixed number of nodes M=400 and level of discretization N=128. In the estimated density plot, we can see that our numerical result for this example is in agreement with the theoretical result in Lindsay and Mallet [34,35], which states that the maximum likelihood estimator can be found in the class of discrete distributions with at most m support points, where Mm. In this example, since we had nine drinking episodes, the estimated density plot displays the support points among 400 nodes for q=(q1,q2).

    In addition, for this sample, the sample mean of q=(q1,q2) is calculated to be

    ˉq=(0.6003,1.2452),

    the sample covariance matrix is calculated to be

    Sq=(0.07060.02640.02640.0483)

    and the sample correlation is calculated to be ρq=0.4519. Based on this, we observe that for our training population consisting of nine drinking episodes, there is a moderate negative association between the parameters q1 and q2. Recall that the parameters q1 and q2 in the system (1.1)–(1.5) represent the normalized diffusivity and the normalized flux gain at the boundary between the dermal and epidermal layers, respectively. It is important to note that the negative correlation is specifically for these nine drinking episodes, and it should not be generalized as the behavior expected in a different sample. In particular, in the next example we have a moderate positive correlation. In addition, it is possible to argue heuristically why one would expect either negative or positive correlation. Significantly more data and testing would be required for a statistical analysis with sufficient power to determine if any conclusion regarding the correlation between the parameters could be drawn.

    We applied the LOOCV method as explained above to the nine drinking episodes from the SCRAM biosensor. Figure 4 shows the measured TAC (i.e., measured by the SCRAM alcohol biosensor) and the estimated TAC (i.e., obtained from our algorithm) for all nine drinking episodes left out in the testing set in the partitioning process, and the conservative 95 % error band for a fixed number of nodes M=400 and level of discretization N=128.

    Figure 4.  The measured TAC, the estimated TAC, and the conservative 95 % error band for nine drinking episodes from the testing set collected using the SCRAM alcohol biosensor using the LOOCV. The figure illustrates a reasonable fit. The second peak in the middle panel of the top row was due to the individual leaving the lab on a warm day, which results in a burst of TAC. This error in the measurement was ignored by the model since this particular test episode was trained based on the remaining eight drinking episodes, which shows that TAC goes to zero as time increases.

    Alcohol researchers and clinicians are particularly interested in certain statistics associated with drinking episodes: The maximum or peak value of the TAC curve, the time at which the peak value of the TAC is attained and the area under the TAC curve. The area under the curve (AUC) is a quantifying measure of exposure to the alcohol that integrates the transdermal alcohol concentration across time. Tables 46 display these statistics along with the measured (or actual) value obtained by the SCRAM alcohol biosensor, as well as the conservative 95 error band for the nine drinking episodes from the testing set. From these tables, we can observe that the 95 error bands do a reasonably good job of capturing the actual values of these statistics, especially for the value of the peak TAC and the area under the curve displayed in Tables 4 and 6. In Figure 4, we can see that there are minor fluctuations in the measured TAC curve. If we smooth the measured TAC curve, we will obtain better results in lowering the error between the time at which the peak value of the TAC is attained using the smoothed TAC curve and the estimated peak time compared to the results obtained in Table 5. The smoothed version results are included in Table A1 in the Appendix.

    Table 4.  The measured peak TAC, estimated peak TAC and the 95 % error band for the nine drinking episodes from the testing set collected using the SCRAM alcohol biosensor.
    Drinking episode Measured peak TAC Estimated peak TAC 95% Error band
    1 0.0585 0.0413 (0.0327, 0.0539)
    2 0.0477 0.0433 (0.0324, 0.0543)
    3 0.0464 0.0391 (0.0320, 0.0501)
    4 0.0417 0.0375 (0.0301, 0.0489)
    5 0.0535 0.0618 (0.0497, 0.0781)
    6 0.0419 0.0361 (0.0287, 0.0465)
    7 0.0450 0.0441 (0.0346, 0.0558)
    8 0.0405 0.0430 (0.0345, 0.0542)
    9 0.0391 0.0402 (0.0313, 0.0508)

     | Show Table
    DownLoad: CSV
    Table 5.  The measured peak time, estimated peak time and the 95 % error band for the nine drinking episodes from the testing set collected using the SCRAM alcohol biosensor.
    Drinking episode Measured peak time Estimated peak time 95% Error band
    1 2.5600 2.4480 (1.9200, 2.8800)
    2 2.2400 2.6080 (2.2400, 2.8800)
    3 2.2400 2.8928 (2.5600, 3.2000)
    4 4.1600 2.9056 (2.5600, 3.2000)
    5 2.2400 2.5728 (2.2400, 2.8800)
    6 2.2400 2.8928 (2.5600, 3.2000)
    7 3.2000 2.9696 (2.5600, 3.2000)
    8 2.8800 2.9312 (2.5600, 3.2000)
    9 3.2000 2.9088 (2.5600, 3.2000)

     | Show Table
    DownLoad: CSV
    Table 6.  The measured AUC, estimated AUC and 95 % error band for nine drinking episodes from the testing set collected using the SCRAM alcohol biosensor.
    Drinking episode Measured AUC Estimated AUC 95% Error band
    1 0.1909 0.1784 (0.1382, 0.2229)
    2 0.1876 0.1799 (0.1343, 0.2321)
    3 0.1868 0.1751 (0.1364, 0.2421)
    4 0.1765 0.1735 (0.1349, 0.2394)
    5 0.2231 0.2963 (0.2203, 0.3911)
    6 0.1151 0.1496 (0.1117, 0.1982)
    7 0.1474 0.1912 (0.1444, 0.2563)
    8 0.1277 0.1898 (0.1412, 0.2505)
    9 0.1493 0.1750 (0.1307, 0.2319)

     | Show Table
    DownLoad: CSV

    For the second example, we applied the LOOCV method to the nine drinking episodes from the WrisTASTM7 alcohol biosensor. Figure 5 shows the measured TAC (i.e., measured by the WrisTASTM7 alcohol biosensor) and the estimated TAC (i.e., obtained from our algorithm) for all nine drinking episodes left out in the testing set in the partitioning process, and the conservative 95 % error band for a fixed number of nodes M=400 and level of discretization N=128. We can observe that we obtained similar results as the previous example which illustrates the consistency of the developed algorithm across different alcohol biosensors.

    Figure 5.  The measured TAC, the estimated TAC and the conservative 95 % error band for nine drinking episodes from the testing set collected using the WrisTASTM7 alcohol biosensor using the LOOCV. The figure illustrates a reasonable fit.

    In fact, Figure 5 shows a better fit than Figure 4, and this might be due to the fact that the WrisTASTM7 biosensor collected data more frequently (i.e., in a smaller time intervals) which resulted in more data points for the training purposes of the model. Another reason may be the difference between the accuracy of the biosensors, which cannot be tested using the dataset used in this paper due to several confounding factors since the two collections of nine drinking episodes in each example do not represent the same drinking episodes; that is, we do not have data collection by the two biosensors under the same conditions to eliminate the confounding factors. If the scientific question that is being addressed is the accuracy of the biosensors, then the researchers should collect data using the two biosensors simultaneously during the same drinking episode to eliminate the confounding factors.

    Similar to the previous example using the SCRAM biosensor, in this example, using the WrisTASTM7 biosensor, we have computed some statistics as the means of comparison between the two datasets. The sample mean of q=(q1,q2) is calculated to be

    ˉq=(0.5696,0.8608),

    the sample covariance matrix is calculated to be

    Sq=(0.02700.00670.00670.0074),

    and the sample correlation is calculated to be ρq=0.4771. We observe a moderate positive association between the parameters q1 and q2, while in the previous example we observed a moderate negative association between the parameters. The statistics obtained in the second example (i.e., using the WrisTASTM7 biosensor) may not be comparable to the statistics obtained in the first example (i.e., using the SCRAM biosensor) due to the fact that the two sets of nine drinking episodes are not identical; in other words, we lack data gathered by the two biosensors in identical conditions to eliminate the potential confounding factors. In addition, the sample size of nine drinking episodes is a small sample size to make any general hypothesis regarding the correlation between the two parameters q1 and q2.

    Tables 79 display these statistics along with the measured (or actual) value obtained by the WrisTASTM7 alcohol biosensor, as well as the conservative 95 error band for the nine drinking episodes from the testing set. From these tables, we can observe that the 95 error bands do a reasonably good job of capturing the actual values of these statistics, especially for the value of the peak TAC and the area under the curve displayed in Tables 7 and 9. In Figure 5, we can see that there are minor fluctuations in the measured TAC curve. Similar to the previous example, if we smooth the measured TAC curve, we will obtain better results in lowering the error between the time at which the peak value of the TAC is attained using the smoothed TAC curve and the estimated peak time compared to the results obtained in Table 8. The smoothed version results are included in Table A2 in the Appendix.

    Table 7.  The measured peak TAC, estimated peak TAC and the 95 % error band for the nine drinking episodes from the testing set collected using the WrisTASTM7 alcohol biosensor.
    Drinking episode Measured peak TAC Estimated peak TAC 95% Error band
    1 0.0228 0.0265 (0.0215, 0.0336)
    2 0.0259 0.0258 (0.0221, 0.0319)
    3 0.0270 0.0218 (0.0195, 0.0266)
    4 0.0222 0.0212 (0.0173, 0.0263)
    5 0.0249 0.0205 (0.0168, 0.0264)
    6 0.0276 0.0263 (0.0221, 0.0337)
    7 0.0287 0.0276 (0.0228, 0.0348)
    8 0.0182 0.0236 (0.0197, 0.0286)
    9 0.0284 0.0296 (0.0247, 0.0371)

     | Show Table
    DownLoad: CSV
    Table 8.  The measured peak time, estimated peak time and the 95 % error band for the nine drinking episodes from the testing set collected using the WrisTASTM7 alcohol biosensor.
    Drinking episode Measured peak time Estimated peak time 95% Error band
    1 2.1667 2.2758 (1.9167, 2.6667)
    2 1.6000 3.0125 (2.8333, 3.5000)
    3 1.6944 2.5542 (2.0833, 3.3333)
    4 2.4167 2.5083 (2.3333, 2.9167)
    5 2.4167 2.3117 (2.1229, 2.5833)
    6 2.4167 2.6917 (2.5000, 3.0000)
    7 1.7500 2.3358 (2.0833, 2.8333)
    8 2.5833 1.9383 (1.8333, 2.0833)
    9 1.5208 2.6233 (2.3333, 3.0833)

     | Show Table
    DownLoad: CSV
    Table 9.  The measured AUC, estimated AUC and 95 % error band for nine drinking episodes from the testing set collected using the WrisTASTM7 alcohol biosensor.
    Drinking episode Measured AUC Estimated AUC 95% Error band
    1 0.0755 0.1054 (0.0912, 0.1229)
    2 0.0937 0.1205 (0.1057, 0.1438)
    3 0.1091 0.1022 (0.0927, 0.1240)
    4 0.0653 0.0876 (0.0759, 0.1023)
    5 0.0902 0.0730 (0.0647, 0.0871)
    6 0.0942 0.1044 (0.0917, 0.1248)
    7 0.1036 0.1213 (0.1063, 0.1448)
    8 0.0716 0.0876 (0.0760, 0.1038)
    9 0.0942 0.1298 (0.1138, 0.1549)

     | Show Table
    DownLoad: CSV

    In this paper, we considered the nonparametric fitting of a population model for the transdermal transport of alcohol based on a maximum likelihood approach applied to a mixed effects statistical model. In estimating a population model, we were actually estimating the distribution of the model parameters and, consequently, the maximum likelihood estimation (MLE) problem was formulated as an optimization problem over a space of feasible probability measures endowed with the weak topology induced by the Prohorov metric. A notable advantage of our developed methodology is that by considering a nonparametric approach, there is no need to make any assumptions about the distribution being estimated.

    By using a first principles physics-based model in the form of a one-dimensional diffusion equation, we were able to capture the essential features of transdermal transport while keeping the dimension of the parameter space low. In this way, we were able to avoid having to introduce regularization so as to mitigate ill-posedness and overfitting. On the other hand, the fact that the model was infinite dimensional being based on a partial differential equation, computing the MLE necessitated finite dimensional approximation.

    We were able to first theoretically demonstrate the existence and then the consistency of our MLE using a decades old result from the literature. The consistency result is with respect to the uncertainty across subjects. It is likely that the consistency results proved in [22] and [21], in the context of a naive pooled statistical model based on a nonlinear least squares estimator, for problems either the same as, or very similar to the one we consider here, would apply for the uncertainty within each subject (i.e., as the resolution of the data with respect to time increases). At present, this is just a hypothesis and a possible avenue for future research; as of yet, we have not carefully examined this possibility.

    In addition, we were able to use linear semigroup theory, in particular the Trotter-Kato theorem, the properties of the weak topology and the Prohorov metric on the space of feasible probability measures to establish a convergence result with respect to the MLE for the finite dimensional approximating estimation problems and the MLE for the estimation problem posed in terms of the original underlying infinite dimensional model.

    We were able to demonstrate the efficacy of our theoretical results numerically first on an example involving simulated data, and then on one involving actual human subject data from an National Institute of Health (NIH) funded study. We used our scheme to obtain the joint density and distribution of the parameters, as well as estimates and conservative 95 error bands for the TAC signal and a number of TAC related statistics of particular interest to researchers and clinicians who work in the area of alcohol use disorder.

    A similar methodology has been developed in [46] using a method called the nnparametric adaptive grid (NPAG) algorithm. In their approach, the grid of support points is iteratively changed in the algorithm, which may add to the computational cost. Also as mentioned in [46], condensing the grid around the nodes with the highest probability may result in estimating only a local maximum likelihood, which may not capture the "true'' distribution to be estimated in the case of a multimodal distribution. One additional difference is that their method was tested on simulated data using a normal distribution for the parameters and it has not been tested on a human subject data. Instead of condensing the grid, we use a fixed grid and directly maximized the log-likelihood function using a mixed-effects model which was tested on both simulated data and on human subject data.

    In addition to consistency with respect to the intrinsic uncertainty, other extensions we are currently looking at include the development of a general framework for estimating random parameters in general finite or infinite dimensional, continuous or discrete-time dynamical systems (e.g., ODEs (ordinary differential equations), PDEs, FDEs (functional differential equations), DEs (difference equations), and more) that would potentially subsume the results presented here as well as in [22] and [21].

    Finally, although correlated, TAC and BrAC are not quantitatively equivalent. BrAC is a good indicator of BAC, which is an informative quantitative measure of the level of alcohol in the body at the moment that is affecting judgment/decision-making, motor coordination, and cognition. On the other hand, TAC is temporally delayed from BrAC since it must pass through the skin prior to being measured, and TAC levels vary across individuals, devices, and within and between drinking episodes, unlike BrAC, which is relatively consistent across people and conditions. This makes TAC difficult to interpret in relation to alcohol levels in the blood at a given time, limiting its usefulness as a quantitative measure of alcohol in the body. Thus, having a way to estimate BrAC from TAC would increase its utility for researchers, clinicians, and individuals who want quantitatively and temporally meaningful information about alcohol levels in the body.

    Since the actual motivation for this investigation is the development of schemes for converting biosensor measured TAC into BAC/BrAC, the next step would be to examine how well population models, estimated using the approach we have presented here, perform when used as part of a scheme that deconvolves an estimate for BAC/BrAC from the TAC signal. In particular, we are interested in comparing it to the schemes, used for this same purpose, developed and implemented in [11,12,28,29,30]. In addition, we are also interested in examining how our uncertainty quantification scheme for the TAC to BAC/BrAC conversion problem performs when compared to the non-physics-based, machine learning inspired schemes developed in [18,19,20,33].

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    We thank the Luczak laboratory students and staff members, particularly Emily Saldich, for their assistance with data collection and management for the SCRAM biosensor. We also thank Dr. Tamara Wall for providing the data for the WrisTASTM7 biosensor.

    This study was funded in part by the National Institute on Alcohol Abuse and Alcoholism (Grant Numbers: R21AA017711 and R01AA026368, S.E.L. and I.G.R.) and by support from the USC Women in Science and Engineering (WiSE) program (L.A.).

    The authors declare no conflicts of interest.

    As mentioned in Section 7.2, smoothing the measured TAC curve results in lowering the error between the time at which the peak value of the TAC is attained using the smoothed TAC curve and the estimated peak time. Here, we include the results using the smoothed version of the TAC curve. The following Tables A1 and A2 correspond to Tables 5 and 8, respectively, with the smoothed TAC curve used instead of the measured TAC curve.

    Table A1.  The peak time using the smoothed TAC curve, estimated peak time, and the 95 error band for the nine drinking episodes from the testing set collected using the SCRAM alcohol biosensor.
    Drinking episode Smoothed peak time Estimated peak time 95% Error band
    1 3.2000 2.4480 (1.9200, 2.8800)
    2 2.8800 2.6080 (2.2400, 2.8800)
    3 2.8800 2.8928 (2.5600, 3.2000)
    4 3.8400 2.9056 (2.5600, 3.2000)
    5 3.2000 2.5728 (2.2400, 2.8800)
    6 2.8800 2.8928 (2.5600, 3.2000)
    7 3.5200 2.9696 (2.5600, 3.2000)
    8 3.5200 2.9312 (2.5600, 3.2000)
    9 3.2000 2.9088 (2.5600, 3.2000)

     | Show Table
    DownLoad: CSV
    Table A2.  The peak time using the smoothed TAC curve, estimated peak time, and the 95 % error band for the nine drinking episodes from the testing set collected using the WrisTASTM7 alcohol biosensor.
    Drinking episode Smoothed peak time Estimated peak time 95% Error band
    1 2.6667 2.2758 (1.9167, 2.6667)
    2 2.1667 3.0125 (2.8333, 3.5000)
    3 2.3333 2.5542 (2.0833, 3.3333)
    4 2.9167 2.5083 (2.3333, 2.9167)
    5 2.6667 2.3117 (2.1229, 2.5833)
    6 2.6667 2.6917 (2.5000, 3.0000)
    7 2.0833 2.3358 (2.0833, 2.8333)
    8 2.6667 1.9383 (1.8333, 2.0833)
    9 2.0833 2.6233 (2.3333, 3.0833)

     | Show Table
    DownLoad: CSV

    To smooth the TAC curve, we used cubic splines and we measured the discrepancy between the smoothed peak time and the estimated peak time using the mean squared error. For the SCRAM biosensor, the mean squared error decreased from 0.314 (prior to smoothing) to 0.293 (after smoothing). For the WrisTASTM7 biosensor, the mean squared error decreased from 0.535 (prior to smoothing) to 0.233 (after smoothing).



    [1] G. Haas, A. Hosmalin, F. Hadida, J. Duntze, P. Debre, B. Autran, Dynamics of HIV variants and specific cytotoxic T-cell recognition in nonprogressors and progressors, Immunol. Lett., 57 (1997), 63–68. https://doi.org/10.1016/S0165-2478(97)00076-X doi: 10.1016/S0165-2478(97)00076-X
    [2] F. Kirchhoff, IV life cycle: Overview. In: T. Hope, M. Stevenson, D. Richman, (eds), Encycl. AIDS, Springer, New York, (2013), 1–9. https://doi.org/10.1007/978-1-4614-9610-6_60-1
    [3] M. A. Nowak, S. Bonhoeffer, G. M. Shaw, R. M. May, Anti-viral drug treatment: Dynamics of resistance in free virus and infected cell populations, J. Theor. Biol., 184 (1997), 203–217. https://doi.org/10.1006/jtbi.1996.0307 doi: 10.1006/jtbi.1996.0307
    [4] T. B. Kepler, A. S. Perelson, Drug concentration heterogeneity facilitates the evolution of drug resistance, Proc. Natl. Acad. Sci. USA, 95 (1998), 11514–11519. https://doi.org/10.1073/pnas.95.20.11514 doi: 10.1073/pnas.95.20.11514
    [5] R. J. Smith, L. M. Wahl, Distinct effects of protease and reverse transcriptase inhibition in an immunological model of HIV-1 infection with impulsive drug effects, Bull. Math. Biol., 66 (2004), 1259–1283. https://doi.org/10.1016/j.bulm.2003.12.004 doi: 10.1016/j.bulm.2003.12.004
    [6] T. H. Zha, O. Castillo, H. Jahanshahi, A. Yusuf, M. O. Alassafi, F. E. Alsaadi, et al., A fuzzy-based strategy to suppress the novel coronavirus (2019-NCOV) massive outbreak, Appl. Comput. Math., 20 2021,160–176.
    [7] M. A. Iqbal, Y. Wang, M. M. Miah, M. S. Osman, Study on DateJimbo-Kashiwara-Miwa equation with conformable derivative dependent on time parameter to find the exact dynamic wave solutions, Fractal Fract., 6 (2022), 1–12. https://doi.org/10.3390/fractalfract6010004 doi: 10.3390/fractalfract6010004
    [8] Y. M. Chu, S. Bashir, M. Ramzan, M. Y. Malik, Model-based comparative study of magnetohydrodynamics unsteady hybrid nanofluid flow between two infinite parallel plates with particle shape effects, Math. Methods Appl. Sci., 2022. http://dx.doi.org/10.1002/mma.8234
    [9] M. Nazeer, F. Hussain, M. I. Khan, A. U. Rehman, E. R. El-Zahar, Y. M. Chu, et al., Theoretical study of MHD electro-osmotically flow of third-grade fluid in micro channel, Appl. Math. Comput., 420 (2022), 126868. https://doi.org/10.1016/j.amc.2021.126868 doi: 10.1016/j.amc.2021.126868
    [10] Y. M. Chu, B. M. Shankaralingappa, B. J. Gireesha, F. Alzahrani, M. I. Khan, S. U. Khan, Combined impact of Cattaneo-Christov double diffusion and radiative heat flux on bio-convective flow of Maxwell liquid configured by a stretched nanomaterial surface, Appl. Math. Comput., 419 (2022), 126883. https://doi.org/10.1016/j.amc.2021.126883 doi: 10.1016/j.amc.2021.126883
    [11] T. H. Zhao, M. I. Khan, Y. M. Chu, Artificial neural networking (ANN) analysis for heat and entropy generation in flow of non-Newtonian fluid between two rotating disks, Math. Methods Appl. Sci., 2021. https://doi.org/10.1002/mma.7310
    [12] L. Wang, M. Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T cells, Math. Biosci., 200 (2006), 44–57. https://doi.org/10.1016/j.mbs.2005.12.026 doi: 10.1016/j.mbs.2005.12.026
    [13] P. K. Srivastava, M. Banerjee, P. Chandra, Modeling the drug therapy for HIV infection, J. Biol. Syst., 17 (2009), 213–223. https://doi.org/10.1142/S0218339009002764 doi: 10.1142/S0218339009002764
    [14] A. Dutta, P. K. Gupta, A mathematical model for transmission dynamics of HIV/AIDS with effect of weak CD4+ T cells, Chin. J. Phys., 56 (2018), 1045–1056. https://doi.org/10.1016/j.cjph.2018.04.004 doi: 10.1016/j.cjph.2018.04.004
    [15] L. Rong, M. A. Gilchrist, Z. Feng, A. S. Perelson, Modeling within host HIV-1 dynamics and the evolution of drug resistance: Trade offs between viral enzyme function and drug susceptibility, J. Theor. Biol., 247 (2007), 804–818. https://doi.org/10.1016/j.jtbi.2007.04.014 doi: 10.1016/j.jtbi.2007.04.014
    [16] Z. Mukandavire, C. Chiyaka, W. Garira, G. Musuka, Mathematical analysis of a sex-structured HIV/AIDS model with a discrete time delay, Nonlinear Anal., Theory Methods Appl., 71 (2009), 1082–1093. https://doi.org/10.1016/j.na.2008.11.026 doi: 10.1016/j.na.2008.11.026
    [17] M. F. Tabassum, M. Saeed, A. Akgül, M. Farman, N. A. Chaudhry, Treatment of HIV/AIDS epidemic model with vertical transmission by using evolutionary Pade-approximation, Chaos Solitons Fractals, 134 (2020), 109686. https://doi.org/10.1016/j.chaos.2020.109686 doi: 10.1016/j.chaos.2020.109686
    [18] A. S. Perelson, Modeling the interaction of the immune system with HIV, In: C. Castillo-Chavez (eds) Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, Springer, Berlin, Heidelberg, (1989), 350–370. https://doi.org/10.1007/978-3-642-93454-4_17
    [19] A. S. Perelson, D. E. Kirschner, R. D. Boer, Dynamics of HIV infection of CD4+ T cells, Math. Biosci., 114 (1993), 81–125. https://doi.org/10.1016/0025-5564(93)90043-A doi: 10.1016/0025-5564(93)90043-A
    [20] S. Arshad, O. Defterli, D. Baleanu, A second order accurate approximation for fractional derivatives with singular and non-singular kernel applied to a HIV model, Appl. Math. Comput., 374 (2020), 125061. https://doi.org/10.1016/j.amc.2020.125061 doi: 10.1016/j.amc.2020.125061
    [21] D. Baleanu, H. Mohammadi, S. Rezapour, Analysis of the model of HIV-1 infection of CD4+ T-cell with a new approach of fractional derivative, Adv. Differ. Equation, 2020 (2020), 1–17. https://doi.org/10.1186/s13662-020-02544-w doi: 10.1186/s13662-020-02544-w
    [22] A. Jan, R. Jan, H. Khan, M. S. Zobaer, R. Shah, Fractional-order dynamics of Rift Valley fever in ruminant host with vaccination, Commun. Math. Biol. Neurosci., 2020 (2020), 1–32. https://doi.org/10.28919/cmbn/5017 doi: 10.28919/cmbn/5017
    [23] R. Jan, M. A. Khan, P. Kumam, P. Thounthong, Modeling the transmission of dengue infection through fractional derivatives, Chaos, Solitons Fractals, 127 (2019), 189–216. https://doi.org/10.1016/j.chaos.2019.07.002 doi: 10.1016/j.chaos.2019.07.002
    [24] F. Wang, M. N. Khan, I. Ahmad, H. Ahmad, H. Abu-Zinadah, Y. M. Chu, Numerical solution of traveling waves in chemical kinetics: time fractional fishers equations, Fractals, 30 (2022), 22400051, 1–11. https://doi.org/10.1142/S0218348X22400515 doi: 10.1142/S0218348X22400515
    [25] Y. M. Chu, A. Ali, M. A. Khan, S. Islam, S. Ullah, Dynamics of fractional order COVID-19 model with a case study of Saudi Arabia, Results Phys., 21 (2021), 103787. https://doi.org/10.1016/j.rinp.2020.103787 doi: 10.1016/j.rinp.2020.103787
    [26] M. A. Khan, S. Ullah, M. Farhan, The dynamics of Zika virus with Caputo fractional derivative, AIMS Math., 4 (2019), 134–146. https://doi.org/10.3934/Math.2019.1.134 doi: 10.3934/Math.2019.1.134
    [27] A. Atangana, Fractal-fractional differentiation and integtion: connecting fractal calculus and fractional calculus to predict complex system, Chaos, Solitons Fractals, 102 (2017), 396–406. https://doi.org/10.1016/j.chaos.2017.04.027 doi: 10.1016/j.chaos.2017.04.027
    [28] S. Bushnaq, S. A. Khan, K. Shah, G. Zaman, Mathematical analysis of HIV/AIDS infection model with Caputo-Fabrizio fractional derivative, Cogent Math. Stat., 5 (2018), 1432521. https://doi.org/10.1080/23311835.2018.1432521 doi: 10.1080/23311835.2018.1432521
    [29] S. Ahmad, A. Ullah, M. Arfan, K. Shah, On analysis of the fractional mathematical model of rotavirus epidemic with the effects of breastfeeding and vaccination under Atangana Baleanu (AB) derivative, Chaos, Solitons Fractals, 140 (2020), 110233. https://doi.org/10.1016/j.chaos.2020.110233 doi: 10.1016/j.chaos.2020.110233
    [30] P. A. Naik, M. Yavuz, J. Zu, The role of prostitution on HIV transmission with memory: A modeling approach, Alexandria Eng. J., 59 (2020), 2513–2531. https://doi.org/10.1016/j.aej.2020.04.016 doi: 10.1016/j.aej.2020.04.016
    [31] Z. Shah, R. Jan, P. Kumam, W. Deebani, M. Shutaywi, Fractional dynamics of HIV with source term for the supply of new CD4+ T-cells depending on the viral load via Caputo-Fabrizio derivative, Molecules, 26 (2021), 1806. https://doi.org/10.3390/molecules26061806 doi: 10.3390/molecules26061806
    [32] J. Kongson, C. Thaiprayoon, W. Sudsutad, Analysis of a fractional model for HIV CD4+ T-cells with treatment under generalized Caputo fractional derivative, AIMS Math., 6 (2021), 7285–7304. https://doi.org/10.3934/math.2021427 doi: 10.3934/math.2021427
    [33] Fatmawati, M. A. Khan, The dynamics of dengue in fection through fractal-fractional operator with real statistical data, Alexandria Eng. J., 60 (2021), 321–336. https://doi.org/10.1016/j.aej.2020.08.018 doi: 10.1016/j.aej.2020.08.018
    [34] P. A. Naik, M. Ghoreishi, J. Zu, Approximate solution of a nonlinear fractional-order HIV model using homotopy analysis method, Int. J. Numer. Anal. Mod., 19 (2022), 52–84.
    [35] D. Baleanu, A. Jajarmi, H. Mohammadi, S. Rezapour, A new study on the mathematical modelling of human liver with Caputo-Fabrizio fractional derivative, Chaos Solitons Fractals, 134 (2020), 109705. https://doi.org/10.1016/j.chaos.2020.109705 doi: 10.1016/j.chaos.2020.109705
    [36] E. Uçara, N. Özdemir, A fractional model of cancer-immune system with Caputo and Caputo-Fabrizio derivatives, Eur. Phys. J. Plus, 136 (2021), 1–17. https://doi.org/10.1140/epjp/s13360-020-00966-9 doi: 10.1140/epjp/s13360-020-00966-9
    [37] Z. Ali, F. Rabiei, K. Shah, T. Khodadadi, Fractal-fractional order dynamical behavior of an HIV/AIDS epidemic mathematical model, Eur. Phys. J. Plus, 136 (2021), 1–17. https://doi.org/10.1140/epjp/s13360-020-00994-5 doi: 10.1140/epjp/s13360-020-00994-5
    [38] S. Ahmad, A. Ullah, A. Akgül, M. D. L. Sen, Study of HIV disease and its association with immune cells under nonsingular and nonlocal fractal-fractional operator, Alexandria Eng. J., (2021), 1904067. https://doi.org/10.1155/2021/1904067
    [39] P. A. Naik, J. Zu, K. M. Owolabi, Modeling the mechanics of viral kinetics under immune control during primary infection of HIV-1 with treatment in fractional order, Phys. A: Stat. Mech. Appl., 545 (2020), 123816. https://doi.org/10.1016/j.physa.2019.123816 doi: 10.1016/j.physa.2019.123816
    [40] P. A. Naik, K. M. Owolabi, M. Yavuz, J. Zu, Chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells, Chaos, Solitons Fractals, 140 (2020), 110272. https://doi.org/10.1016/j.chaos.2020.110272 doi: 10.1016/j.chaos.2020.110272
    [41] P. A. Naik, J. Zu, M. Ghoreishi, Estimating the approximate analytical solution of HIV viral dynamic model by using homotopy analysis method, Chaos, Solitons Fractals, 131 (2020), 109500. https://doi.org/10.1016/j.chaos.2019.109500 doi: 10.1016/j.chaos.2019.109500
    [42] M. M. El-Dessoky, M. A. Khan, Modeling and analysis of an epidemic model with fractal-fractional Atangana-Baleanu derivative, Alexandria Eng. J., 61 (2022), 729–746. https://doi.org/10.1016/j.aej.2021.04.103
    [43] X. Q. Zhao, The theory of basic reproduction ratios, in Dynamical Systems in Population Biology, CMS Books in Mathematics, Springer, Cham, (2017), 285–315. https://doi.org/10.1007/978-3-319-56433-3_11
    [44] C. Castillo-Chavez, Z. Feng, W. Huang, On the computation of R0 and its role on global stability, in Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, Springer, Berlin, 2002. https://doi.org/10.1007/978-1-4757-3667-0_13
    [45] E. Ahmeda, A. S. Elgazzar, On fractional order differential equations model for nonlocal epidemics, Phys. A: Stat. Mech. Appl., 379 (2007), 607–614. https://doi.org/10.1016/j.physa.2007.01.010 doi: 10.1016/j.physa.2007.01.010
    [46] S. Bourafaa, M. S. Abdelouahaba, A. Moussaoui, On some extended Routh–Hurwitz conditions for fractional-order autonomous systems of order α(0,2) and their applications to some population dynamic models, Chaos, Solitons Fractals, 133 (2020), 109623. https://doi.org/10.1016/j.chaos.2020.109623 doi: 10.1016/j.chaos.2020.109623
    [47] A. Granas, J. Dugundji, Fixed Point Theory, Springer, New York, 2003. https://doi.org/10.1007/978-0-387-21593-8
    [48] D. H. Griffel, Applied Functional Analysis, Ellis Horwood: Chichester, UK, 1981.
    [49] A. Atangana, S. I. Araz, New Numerical Scheme with Newton Polynomial: Theory, Methods, and Applications, 1st edition, Elsevier, 2021. https://doi.org/10.1016/C2020-0-02711-8
    [50] A. T. Haase, K. Henry, M. Zupancic, G. Sedgewick, R. A. Faust, H. Melroe, et al., Quantitative image analysis of HIV-1 infection in lymphoid tissue, Science, 274 (1996), 985–989. https://doi.org/10.1126/science.274.5289.9 doi: 10.1126/science.274.5289.9
    [51] R. D. Hockett, J. M. Kilby, C. A. Derdeyn, M. S. Saag, M. Sillers, K. Squires, et al., Constant mean viral copy number per infected cell in tissues regardless of high, low, or undetectable plasma HIV RNA, J. Exp. Med., 189 (1999), 1545–1554. https://doi.org/10.1084/jem.189.10.1545 doi: 10.1084/jem.189.10.1545
    [52] P. Essunger, A. S. Perelson, Modeling HIV infection of CD4+ T-cell subpopulations, J. Theor. Biol., 170 (1994), 367–391. https://doi.org/10.1006/jtbi.1994.1199 doi: 10.1006/jtbi.1994.1199
    [53] H. Mohri, S. Bonhoeffer, S. Monard, A. S. Perelson, D. D. Ho, Rapid turnover of T lymphocytes in SIV-infected rhesus macaques, Science, 279 (1998), 1223–1227. https://doi.org/10.1126/science.279.5354.1223 doi: 10.1126/science.279.5354.1223
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2259) PDF downloads(126) Cited by(5)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog