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

Computational study of a homogenized nonlinear generalization of Timoshenko beam proposed by Turco et al.

  • To my mentor and friend Emilio Turco, in the occasion of his 60th birthday.
    EB
  • Mechanical metamaterials are most often assemblies of stocky beam elements connected through rigid connections, hinges, or flexural joints. The description of these materials through classical beam theories is challenging because of the wide variety of complex phenomena observed in the severe deformation regime mechanical metamaterials must undergo and because most classical beam theories can only be applied to elements with sufficiently high slenderness. In the spirit of Hencky, Turco et al. (2020) has recently formulated an intrinsically discrete nonlinear elastic model suitable for the design of mechanical metamaterials. The objective of this contribution was to present a numerical study of the nonlinear generalization of the Timoshenko beam that results from the asymptotic homogenization of the discrete model introduced by Turco et al. The present numerical study took into account several loading cases and elucidated the sensitivity of the homogenized continuum with respect to axial, bending, and shear stiffness parameters, as well as to load imperfections, in terms of mechanical behavior, including buckling onset and post-critical behavior. It was found that the predictions obtained with the homogenized model in the large deformation regime matched excellently with those of the discrete model proposed by Turco et al.

    Citation: Jose Manuel Torres Espino, Emilio Barchiesi. Computational study of a homogenized nonlinear generalization of Timoshenko beam proposed by Turco et al.[J]. Networks and Heterogeneous Media, 2024, 19(3): 1133-1155. doi: 10.3934/nhm.2024050

    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
  • Mechanical metamaterials are most often assemblies of stocky beam elements connected through rigid connections, hinges, or flexural joints. The description of these materials through classical beam theories is challenging because of the wide variety of complex phenomena observed in the severe deformation regime mechanical metamaterials must undergo and because most classical beam theories can only be applied to elements with sufficiently high slenderness. In the spirit of Hencky, Turco et al. (2020) has recently formulated an intrinsically discrete nonlinear elastic model suitable for the design of mechanical metamaterials. The objective of this contribution was to present a numerical study of the nonlinear generalization of the Timoshenko beam that results from the asymptotic homogenization of the discrete model introduced by Turco et al. The present numerical study took into account several loading cases and elucidated the sensitivity of the homogenized continuum with respect to axial, bending, and shear stiffness parameters, as well as to load imperfections, in terms of mechanical behavior, including buckling onset and post-critical behavior. It was found that the predictions obtained with the homogenized model in the large deformation regime matched excellently with those of the discrete model proposed by Turco et al.



    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] J. Alibert, P. Seppecher, F. dell'Isola, Truss modular beams with deformation energy depending on higher displacement gradients, Math Mech Solids, 8 (2023), 51–73. https://doi.org/10.1177/1081286503008001658 doi: 10.1177/1081286503008001658
    [2] H. Altenbach, V. Eremeyev, On the constitutive equations of viscoelastic micropolar plates and shells of differential type, Math. Mech. Complex Syst., 3 (2015), 273–283. https://doi.org/10.2140/memocs.2015.3.273 doi: 10.2140/memocs.2015.3.273
    [3] E. Barchiesi, Multi-scale and multi-physics: towards next-generation engineering materials, Continuum Mech. Thermodyn., 32 (2020), 541–554.
    [4] E. Barchiesi, F. dell'Isola, A. M. Bersani, E. Turco, Equilibria determination of elastic articulated duoskelion beams in 2d via a riks-type algorithm, Int. J. Non Linear Mech, 128 (2021), 103628. https://doi.org/10.1016/j.ijnonlinmec.2020.103628 doi: 10.1016/j.ijnonlinmec.2020.103628
    [5] E. Barchiesi, F. dell'Isola, F. Hild, On the validation of homogenized modeling for bi-pantographic metamaterials via digital image correlation, Int J Solids Struct, 208 (2021), 49–62. https://doi.org/10.1016/j.ijsolstr.2020.09.036 doi: 10.1016/j.ijsolstr.2020.09.036
    [6] K. Barri, Q. Zhang, J. Kline, W. Lu, J. Luo, Z. Sun, et al., Multifunctional nanogenerator-integrated metamaterial concrete systems for smart civil infrastructure, Adv. Mater., 35 (2023), 2211027. https://doi.org/10.1002/adma.202211027 doi: 10.1002/adma.202211027
    [7] R. P. Bohara, S. Linforth, T. Nguyen, A. Ghazlan, T. Ngo, Anti-blast and -impact performances of auxetic structures: A review of structures, materials, methods, and fabrications, Eng Struct, 276 (2023), 115377. https://doi.org/10.1016/j.engstruct.2022.115377 doi: 10.1016/j.engstruct.2022.115377
    [8] C. Boutin, F. dell'Isola, Green's functions and integral representation of generalized continua: the case of orthogonal pantographic lattices, Z. Angew. Math. Phys., 72 (2021), 58. https://doi.org/10.1007/s00033-021-01480-3 doi: 10.1007/s00033-021-01480-3
    [9] M. S. Chaki, V. A. Eremeyev, A. K. Singh, Surface and interfacial anti-plane waves in micropolar solids with surface energy, Math Mech Solids, 26 (2021), 708–721. https://doi.org/10.1177/1081286520965646 doi: 10.1177/1081286520965646
    [10] F. Cornacchia, F. Fabbrocino, N. Fantuzzi, R. Luciano, R. Penna, Analytical solution of cross-and angle-ply nano plates with strain gradient theory for linear vibrations and buckling, Mech Adv. Mater. Struct, 28 (2021), 1201–1215. https://doi.org/10.1080/15376494.2019.1655613 doi: 10.1080/15376494.2019.1655613
    [11] H. Darban, R. Luciano, A. Caporale, F. Fabbrocino., Higher modes of buckling in shear deformable nanobeams, Int J Eng Sci, 154 (2020), 103338. https://doi.org/10.1016/j.ijengsci.2020.103338 doi: 10.1016/j.ijengsci.2020.103338
    [12] M. De Angelo, L. Placidi, N. Nejadsadeghi, A. Misra, Non-standard timoshenko beam model for chiral metamaterial: Identification of stiffness parameters, Mech Res Commun, 103 (2020), 103462. https://doi.org/10.1016/j.ijengsci.2020.103338 doi: 10.1016/j.ijengsci.2020.103338
    [13] I. Elishakoff, Who developed the so-called timoshenko beam theory?, Math Mech Solids, 25 (2020), 97–116. https://doi.org/10.1177/1081286519856931 doi: 10.1177/1081286519856931
    [14] V. Eremeyev, Wojciech Pietraszkiewicz, Refined theories of plates and shells, ZAMM-Z Angew Math Me, 94 (2014), 5.
    [15] V. A. Eremeyev, Two-and three-dimensional elastic networks with rigid junctions: modeling within the theory of micropolar shells and solids, Acta Mech, 230 (2019), 3875–3887. https://doi.org/10.1007/s00707-019-02527-3 doi: 10.1007/s00707-019-02527-3
    [16] V. A. Eremeyev, L. P. Lebedev, V. Konopińska-Zmysłowska, On solvability of initial boundary-value problems of micropolar elastic shells with rigid inclusions, Math Mech Solids, 27 (2022), 1800–1812. https://doi.org/10.1177/10812865211073149 doi: 10.1177/10812865211073149
    [17] V. A. Eremeyev, W. Pietraszkiewicz, Material symmetry group and constitutive equations of micropolar anisotropic elastic solids, Math Mech Solids, 21 (2016), 210–221. https://doi.org/10.1177/1081286515582862 doi: 10.1177/1081286515582862
    [18] V. A. Eremeyev, E. Turco, Enriched buckling for beam-lattice metamaterials, Mech Res Commun, 103 (2020), 103458. https://doi.org/10.1016/j.mechrescom.2019.103458 doi: 10.1016/j.mechrescom.2019.103458
    [19] F. Fabbrocino, M. Funari, F. Greco, P. Lonetti, R. Luciano, R. Penna, Dynamic crack growth based on moving mesh method, Compos Part B-eng, 174 (2019), 107053. https://doi.org/10.1016/j.compositesb.2019.107053 doi: 10.1016/j.compositesb.2019.107053
    [20] N. Feng, Y. Tie, S. Wang, J. Guo, A novel 3D bidirectional auxetic metamaterial with lantern-shape: Elasticity aspects and potential for load-bearing structure, Compos Struct, 321 (2023), 117221. https://doi.org/10.1016/j.compstruct.2023.117221 doi: 10.1016/j.compstruct.2023.117221
    [21] M. F. Funari, S. Spadea, F. Fabbrocino, R. Luciano, A moving interface finite element formulation to predict dynamic edge debonding in frp-strengthened concrete beams in service conditions, Fibers, 8 (2020), 42. https://doi.org/10.3390/fib8060042 doi: 10.3390/fib8060042
    [22] I. Giorgio, F. Hild, E. Gerami, F. dell'Isola, A. Misra, Experimental verification of 2D cosserat chirality with stretch-micro-rotation coupling in orthotropic metamaterials with granular motif, Mech Res Commun, 126 (2022), 104020. https://doi.org/10.1016/j.mechrescom.2022.104020 doi: 10.1016/j.mechrescom.2022.104020
    [23] I. Giorgio, V. Varano, F. dell'Isola, N. L. Rizzi, Two layers pantographs: A 2D continuum model accounting for the beams' offset and relative rotations as averages in SO(3) lie groups, Int J Solids Struct, 216 (2021), 43–58. https://doi.org/10.1016/j.ijsolstr.2021.01.018 doi: 10.1016/j.ijsolstr.2021.01.018
    [24] D. Han, X. Ren, Y. Zhang, X. Y. Zhang, X. G. Zhang, C. Luo, Y. M. Xie, Lightweight auxetic tamaterials: Design and characteristic study, Compos Struct, 293 (2022), 115706. https://doi.org/10.1016/j.compstruct.2022.115706 doi: 10.1016/j.compstruct.2022.115706
    [25] H. Hencky, Über die angenäherte Lösung von Stabilitätsproblemen im Raum mittels der elastischen Gelenkkette, (Germany), Doctoral Thesis of W. Engelmann, Leipzig, 1921.
    [26] F. Hild, A. Misra, F. dell'Isola, Multiscale DIC applied to pantographic structures, Exp Mech, 61 (2021), 431–443. https://doi.org/10.1007/s11340-020-00636-y doi: 10.1007/s11340-020-00636-y
    [27] S. Jin, Y. P. Korkolis, Y. Li, Shear resistance of an auxetic chiral mechanical metamaterial, Int J Solid Struct, 174 (2019), 28–37. https://doi.org/10.1016/j.ijsolstr.2019.06.005 doi: 10.1016/j.ijsolstr.2019.06.005
    [28] N. Karathanasopoulos, F. Dos Reis, M. Diamantopoulou, J. F. Ganghoffer, Mechanics of beams made from chiral metamaterials: Tuning deflections through normal-shear strain couplings, Mater Des, 189 (2020), 108520.
    [29] M. Khezri, K. Rasmussen, Functionalising buckling for structural morphing in kinetic façades: Concepts, strategies and applications, Thin Wall Struct, 180 (2022), 109749. https://doi.org/10.1016/j.matdes.2020.108520 doi: 10.1016/j.matdes.2020.108520
    [30] J. A. Kraus, J. J. Rimoli, An elastica theory for compressible imperfect beams with application to mechanical metamaterials, Mech Res Commun, 131 (2023), 104147. https://doi.org/10.1016/j.mechrescom.2023.104147 doi: 10.1016/j.mechrescom.2023.104147
    [31] R. Luciano, H. Darban, C. Bartolomeo, F. Fabbrocino, D. Scorza, Free flexural vibrations of nanobeams with non-classical boundary conditions using stress-driven nonlocal model, Mech Res Commun, 107 (2020), 103536. https://doi.org/10.1016/j.mechrescom.2020.103536 doi: 10.1016/j.mechrescom.2020.103536
    [32] G. Mancusi, F. Fabbrocino, L. Feo, F. Fraternali, Size effect and dynamic properties of 2D lattice materials, Compos Part B-eng, 112 (2017), 235–242. https://doi.org/10.1016/j.compositesb.2016.12.026 doi: 10.1016/j.compositesb.2016.12.026
    [33] P. K. Masjedi, P. M. Weaver, Analytical solution for arbitrary large deflection of geometrically exact beams using the homotopy analysis method, Appl Math Model, 103 (2022), 516–542. https://doi.org/10.1016/j.apm.2021.10.037 doi: 10.1016/j.apm.2021.10.037
    [34] A. Nuñez-Labielle, J. Cante, A. Huespe, J. Oliver, Towards shock absorbing hyperelastic metamaterial design. (i) macroscopic scale: Computational shock-capturing, Comput. Methods Appl. Mech. Eng., 393 (2022), 114732. https://doi.org/10.1016/j.cma.2022.114732 doi: 10.1016/j.cma.2022.114732
    [35] L. Paleari, M. Bragaglia, F. Fabbrocino, R. Luciano, F. Nanni, Self-monitoring performance of 3d-printed poly-ether-ether-ketone carbon nanotube composites, Polymers, 15 (2022), 8. https://doi.org/10.3390/polym15010008 doi: 10.3390/polym15010008
    [36] D. Qi, P. Zhang, W. Wu, K. Xin, H. Liao, Y. Li, et al., Innovative 3D chiral metamaterials under large deformation: Theoretical and experimental analysis, Int J Solids Struct, 202 (2020), 787–797. https://doi.org/10.1016/j.ijsolstr.2020.06.047 doi: 10.1016/j.ijsolstr.2020.06.047
    [37] A. Rafsanjani, K. Bertoldi, A. R. Studart, Programming soft robots with flexible mechanical metamaterials, Sci. Robot., 4 (2019), eaav7874. https://doi.org/10.1126/scirobotics.aav7874 doi: 10.1126/scirobotics.aav7874
    [38] H. Reda, S. Alavi, M. Nasimsobhan, J. Ganghoffer, Homogenization towards chiral cosserat continua and applications to enhanced timoshenko beam theories, Mech Mater, 155 (2021), 103728. https://doi.org/10.1016/j.mechmat.2020.103728 doi: 10.1016/j.mechmat.2020.103728
    [39] E. Riks, The Application of Newton's Method to the Problem of Elastic Stability, J. Appl. Mech., 39 (1972), 1060–1065. https://doi.org/10.1115/1.3422829 doi: 10.1115/1.3422829
    [40] P. Sengsri, S. Kaewunruen, Additive manufacturing meta-functional composites for engineered bridge bearings: A review, Constr Build Mater, 262 (2020), 120535. https://doi.org/10.1115/1.3422829 doi: 10.1115/1.3422829
    [41] A. Skrzat, V. A. Eremeyev, On the effective properties of foams in the framework of the couple stress theory, Continuum Mech. Thermodyn., 32 (2020), 1779–1801. https://doi.org/10.1007/s00161-020-00880-6 doi: 10.1007/s00161-020-00880-6
    [42] J. U. Surjadi, L. Gao, H. Du, X. Li, X. Xiong, N. X. Fang, et al., Mechanical metamaterials and their engineering applications, Adv Eng Mater, 21 (2019), 1800864. https://doi.org/10.1002/adem.201800864 doi: 10.1002/adem.201800864
    [43] E. Turco, E. Barchiesi, I. Giorgio, F. dell'Isola, A lagrangian hencky-type non-linear model suitable for metamaterials design of shearable and extensible slender deformable bodies alternative to timoshenko theory, Int. J. Non Linear Mech, 123 (2020), 103481. https://doi.org/10.1016/j.ijnonlinmec.2020.103481 doi: 10.1016/j.ijnonlinmec.2020.103481
    [44] E. N. Vilchevskaya, W. Müller, V. A. Eremeyev, Extended micropolar approach within the framework of 3M theories and variations thereof, Continuum Mech. Thermodyn., 34 (2022), 533–554. https://doi.org/10.1007/s00161-021-01072-6 doi: 10.1007/s00161-021-01072-6
    [45] J. Wang, S. Zhu, L. Chen, T. Liu, H. Liu, Z. Lv, et al., Data mining from a hierarchical dataset for mechanical metamaterials composed of curved-sides triangles, Compos Struct, 319 (2023), 117153.
    [46] O. Weeger, Numerical homogenization of second gradient, linear elastic constitutive models for cubic 3D beam-lattice metamaterials, Int J Solids Struct, 224 (2021), 111037. https://doi.org/10.1016/j.ijsolstr.2021.03.024 doi: 10.1016/j.ijsolstr.2021.03.024
    [47] W. Wu, P. Liu, Z. Kang, A novel mechanical metamaterial with simultaneous stretching- and compression-expanding property, Mater Des, 208 (2021), 109930. https://doi.org/10.1016/j.matdes.2021.109930 doi: 10.1016/j.matdes.2021.109930
    [48] B. Yang, M. Bacciocchi, N. Fantuzzi, R. Luciano, F. Fabbrocino, Wave propagation in periodic nano structures through second strain gradient elasticity, Int. J. Mech. Sci. Struct, 260 (2023), 108639. https://doi.org/10.1016/j.ijmecsci.2023.108639 doi: 10.1016/j.ijmecsci.2023.108639
    [49] B. Yang, M. Bacciocchi, N. Fantuzzi, R. Luciano, F. Fabbrocino, Computational simulation and acoustic analysis of two-dimensional nano-waveguides considering second strain gradient effects, Comput Struct, 296 (2024), 107299. https://doi.org/10.1016/j.compstruc.2024.107299 doi: 10.1016/j.compstruc.2024.107299
    [50] P. Zhang, D. Qi, R. Xue, K. Liu, W. Wu, Y. Li, Mechanical design and energy absorption performances of rational gradient lattice metamaterials, Compos Struct, 277 (2021), 114606. https://doi.org/10.1016/j.compstruct.2021.114606 doi: 10.1016/j.compstruct.2021.114606
    [51] X. Y. Zhang, X. Ren, Y. Zhang, Y. M. Xie, A novel auxetic metamaterial with enhanced mechanical properties and tunable auxeticity, Thin Wall Struct, 174 (2022), 109162. https://doi.org/10.1016/j.tws.2022.109162 doi: 10.1016/j.tws.2022.109162
    [52] S. Zolfaghari, D. Mostofinejad, N. Fantuzzi, R. Luciano, F. Fabbrocino, Experimental evaluation of FRP-concrete bond using externally-bonded reinforcement on grooves (EBROG) method, Compos Struct, 310 (2023), 116693. https://doi.org/10.1016/j.compstruct.2023.116693 doi: 10.1016/j.compstruct.2023.116693
  • Reader Comments
  • © 2024 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(1055) PDF downloads(47) Cited by(0)

Figures and Tables

Figures(15)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog