Processing math: 100%
Research article Special Issues

Impact of antibody-level on viral shedding in B.1.617.2 (Delta) variant-infected patients analyzed using a joint model of longitudinal and time-to-event data

  • These authors have contributed equally to this work
  • Knowledge of viral shedding remains limited. Repeated measurement data have been rarely used to explore the influencing factors. In this study, a joint model was developed to explore and validate the factors influencing the duration of viral shedding based on longitudinal data and survival data. We divided 361 patients infected with Delta variant hospitalized in Nanjing Second Hospital into two groups (≤ 21 days group and > 21 days group) according to the duration of viral shedding, and compared their baseline characteristics. Correlation analysis was performed to identify the factors influencing the duration of viral shedding. Further, a joint model was established based on longitudinal data and survival data, and the Markov chain Monte Carlo algorithm was used to explain the influencing factors. In correlation analysis, patients having received vaccination had a higher antibody level at admission than unvaccinated patients, and with the increase of antibody level, the duration of viral shedding shortened. The linear mixed-effects model showed the longitudinal variation of logSARS-COV-2 IgM sample/cutoff (S/CO) values, with a parameter estimate of 0.193 and a standard error of 0.017. Considering gender as an influencing factor, the parameter estimate of the Cox model and their standard error were 0.205 and 0.1093 (P = 0.608), the corresponding OR value was 1.228. The joint model output showed that SARS-COV-2 IgM (S/CO) level was strongly associated with the risk of a composite event at the 95% confidence level, and a doubling of SARS-COV-2 IgM (S/CO) level was associated with a 1.38-fold (95% CI: [1.16, 1.72]) increase in the risk of viral non-shedding. A higher antibody level in vaccinated patients, as well as the presence of IgM antibodies in serum, can accelerate shedding of the mutant virus. This study provides some evidence support for vaccine prevention and control of COVID-19 variants.

    Citation: Yi Yin, Ting Zeng, Miao Lai, Zemin Luan, Kai Wang, Yuhang Ma, Zhiliang Hu, Kai Wang, Zhihang Peng. Impact of antibody-level on viral shedding in B.1.617.2 (Delta) variant-infected patients analyzed using a joint model of longitudinal and time-to-event data[J]. Mathematical Biosciences and Engineering, 2023, 20(5): 8875-8891. doi: 10.3934/mbe.2023390

    Related Papers:

    [1] Alessandra Jannelli, Maria Paola Speciale . On the numerical solutions of coupled nonlinear time-fractional reaction-diffusion equations. AIMS Mathematics, 2021, 6(8): 9109-9125. doi: 10.3934/math.2021529
    [2] Sara S. Alzaid, Pawan Kumar Shaw, Sunil Kumar . A study of Ralston's cubic convergence with the application of population growth model. AIMS Mathematics, 2022, 7(6): 11320-11344. doi: 10.3934/math.2022632
    [3] Xiaopeng Yi, Chongyang Liu, Huey Tyng Cheong, Kok Lay Teo, Song Wang . A third-order numerical method for solving fractional ordinary differential equations. AIMS Mathematics, 2024, 9(8): 21125-21143. doi: 10.3934/math.20241026
    [4] Najat Almutairi, Sayed Saber . Chaos control and numerical solution of time-varying fractional Newton-Leipnik system using fractional Atangana-Baleanu derivatives. AIMS Mathematics, 2023, 8(11): 25863-25887. doi: 10.3934/math.20231319
    [5] Sara S. Alzaid, Pawan Kumar Shaw, Sunil Kumar . A numerical study of fractional population growth and nuclear decay model. AIMS Mathematics, 2022, 7(6): 11417-11442. doi: 10.3934/math.2022637
    [6] Ali Turab, Hozan Hilmi, Juan L.G. Guirao, Shabaz Jalil, Nejmeddine Chorfi, Pshtiwan Othman Mohammed . The Rishi Transform method for solving multi-high order fractional differential equations with constant coefficients. AIMS Mathematics, 2024, 9(2): 3798-3809. doi: 10.3934/math.2024187
    [7] Din Prathumwan, Thipsuda Khonwai, Narisara Phoochalong, Inthira Chaiya, Kamonchat Trachoo . An improved approximate method for solving two-dimensional time-fractional-order Black-Scholes model: a finite difference approach. AIMS Mathematics, 2024, 9(7): 17205-17233. doi: 10.3934/math.2024836
    [8] Naimi Abdellouahab, Keltum Bouhali, Loay Alkhalifa, Khaled Zennir . Existence and stability analysis of a problem of the Caputo fractional derivative with mixed conditions. AIMS Mathematics, 2025, 10(3): 6805-6826. doi: 10.3934/math.2025312
    [9] Abdon Atangana, Seda İğret Araz . A successive midpoint method for nonlinear differential equations with classical and Caputo-Fabrizio derivatives. AIMS Mathematics, 2023, 8(11): 27309-27327. doi: 10.3934/math.20231397
    [10] Khalid K. Ali, Mohamed A. Abd El Salam, Mohamed S. Mohamed . Chebyshev fifth-kind series approximation for generalized space fractional partial differential equations. AIMS Mathematics, 2022, 7(5): 7759-7780. doi: 10.3934/math.2022436
  • Knowledge of viral shedding remains limited. Repeated measurement data have been rarely used to explore the influencing factors. In this study, a joint model was developed to explore and validate the factors influencing the duration of viral shedding based on longitudinal data and survival data. We divided 361 patients infected with Delta variant hospitalized in Nanjing Second Hospital into two groups (≤ 21 days group and > 21 days group) according to the duration of viral shedding, and compared their baseline characteristics. Correlation analysis was performed to identify the factors influencing the duration of viral shedding. Further, a joint model was established based on longitudinal data and survival data, and the Markov chain Monte Carlo algorithm was used to explain the influencing factors. In correlation analysis, patients having received vaccination had a higher antibody level at admission than unvaccinated patients, and with the increase of antibody level, the duration of viral shedding shortened. The linear mixed-effects model showed the longitudinal variation of logSARS-COV-2 IgM sample/cutoff (S/CO) values, with a parameter estimate of 0.193 and a standard error of 0.017. Considering gender as an influencing factor, the parameter estimate of the Cox model and their standard error were 0.205 and 0.1093 (P = 0.608), the corresponding OR value was 1.228. The joint model output showed that SARS-COV-2 IgM (S/CO) level was strongly associated with the risk of a composite event at the 95% confidence level, and a doubling of SARS-COV-2 IgM (S/CO) level was associated with a 1.38-fold (95% CI: [1.16, 1.72]) increase in the risk of viral non-shedding. A higher antibody level in vaccinated patients, as well as the presence of IgM antibodies in serum, can accelerate shedding of the mutant virus. This study provides some evidence support for vaccine prevention and control of COVID-19 variants.



    Fractional calculus is a branch of study that of any order integral or derivative. It is obtained by replacing the integer order of differential equations with the fractional order. Because fractional derivative has heredity or memory, the fractional differential equation is widely used in various fields. In the past few decades, fractional derivative has been studied by a large number of scholars. They firstly proved a (3θ)-order convergence and stability scheme for time-fractional derivative θ under the time step k2 in [1]. The similar (3θ)-order scheme to efficiently solve the time-fractional diffusion equation in time was first constructed in [2]. Based on the idea of block-by-block approach method, the reference [3] presented a general technique to construct high order schemes for the numerical solution of the fractional ordinary differential equations. The reference [4] established a fractional Gronwall inequality to prove some stability and convergence estimates of schemes for fractional reaction-subdiffusion problems. In [5], a space-time finite element method was used to solve the distributed-order time fractional reaction diffusion equations. The scheme with (3θ) convergence order for time step k2 and (2θ) convergence order for time step k=1 was constructed for time-fractional derivative in [6]. They proved the (2θ)-order convergence and stability of the time scheme, where θ is the order of the time-fractional derivative in [7]. A series of high order numerical approximations to θ-Caputo derivatives (0<θ<2) and Riemann-Liouville derivatives were given in [8]. They investigated numerical solutions of coupled nonlinear time-fractional reaction-diffusion equations obtained by the Lie symmetry analysis in [9]. They built a robust finite difference scheme for a nonlinear Caputo fractional differential equation on an adaptive grid with first-order convergent result in [10]. Spectral methods was popularly used to solve fractional derivative in [11,12,13]. An implicit numerical method was constructed to solve the two-dimensional θ(1<θ<2) time fractional diffusion-wave equation with the convergence order 3θ in [14]. They used the piecewise linear and quadratic Lagrange interpolation functions to construct a (3θ)-order numerical approximate method for the Caputo fractional derivative in [15]. In [16], a fast high order numerical method for nonlinear fractional-order differential equations with non-singular kernel was constructed. The numerical scheme for nonlinear fractional differential equations was constructed by the Galerkin finite element method in [17]. The reference [18] established a hp-version Chebyshev collocation method for nonlinear fractional differential equations. The Spectral collocation method for nonlinear Riemann-Liouville fractional differential equations was given in [19]. The above high order schemes are either implicit or low-order schemes are used in the first step or obtained by discretizing the equivalent form of the original equation. Therefore, we will give a high order numerical scheme with uniform convergence order in this paper using the direct numerical discretization of original equations based on the idea of [1,2,3].

    This paper is arranged as follows: In Section 2 the higher order numerical scheme is proposed. And the local truncation error of the constructed higher order numerical scheme is given in Section 3. The convergence analysis is studied in Section 4. In Section 5 some numerical examples are given. Finally, some conclusion are given in Section 6.

    We consider the following initial value problem of nonlinear fractional ordinary differential equations

    {0Dθtu(t)=f(t,u(t)),0tT,0<θ<1,(2.1)u(0)=u0,(2.2)

    where θ is the order of the fractional derivative, 0Dθtu(t) in (2.1) is defined as the Caputo fractional derivatives of order θ given by

    0Dθtu(t)=t0ω1θ  (tτ)u(τ)dτ, (2.3)

    with ω1θ   is defined by ω1θ  (t)=1tθΓ(1θ), and Γ() denotes Gamma function [20].

    First, the finite difference scheme is introduced to discretize the fractional derivative. Dividing the interval [0,T] into K equal subintervals, set tk=kΔt,k=0,1,,K, with Δt=TK, the numerical solution of (2.1) at tk is uk, and 0Dθtu(tk)=0Dθtuk, fk=f(tk,uk).

    Next, we start discretizing the fractional derivative (2.1). Firstly, the values of u(t) on t1 and t2 are determined. Using Quadratic Lagrange interpolation, the approximation formula of u(t) on the interval [t0, t2] is

    J[t0,t2]u(t)=ϖ0,0(t)u0+ϖ1,0(t)u1+ϖ2,0(t)u2, (2.4)

    where ϖi,0(t),i=0,1,2, are the Quadratic Lagrangian interpolation basis functions on point t0, t1 and t2, are defined as following

    ϖ0,0(t)=(tt1)(tt2)2Δt2, ϖ1,0(t)=(tt0)(tt2)Δt2, ϖ2,0(t)=(tt0)(tt1)2Δt2.

    Let Δuk=ukuk1,Δu0=u0,k1. when k=1,2, we use 0Dθt(J[t0,t2]u)(t1), 0Dθt(J[t0,t2]u)(t2) to approximate 0Dθtu(t1), 0Dθtu(t2), respectively, then we get

    0Dθtu(t1)=t10ω1θ  (t1τ)u(τ)dτt10ω1θ  (t1τ)[J[t0,t2]u(τ)]dτ=1Δt2t10ω1θ  (t1τ)[(τt12)Δu2(τt32)Δu1]dτ=B11Δu1+B21Δu2 0DθΔtu1, (2.5)
    0Dθtu(t2)=t20ω1θ  (t2τ)u(τ)dτt20ω1θ  (t2τ)[J[t0,t2]u(τ)]dτ=1Δt2t20ω1θ  (t2τ)[(τt12)Δu2(τt32)Δu1]dτ=B12Δu1+B22Δu2 0DθΔtu2, (2.6)

    where

    B1j=1Δt2tj0ω1θ  (tjτ)(τt32)dτ,  B2j=1Δt2tj0ω1θ  (tjτ)(τt12)dτ,j=1,2. (2.7)

    When k3, we have

    0Dθtu(tk)=t10ω1θ  (tkτ)u(τ)dτ+k1j=1tj+1tjω1θ  (tkτ)u(τ)dτt10ω1θ  (tkτ)[J[t0,t2]u(τ)]dτ+k1j=1tj+1tjω1θ  (tkτ)[J[tj1,tj+1]u(τ)]dτ=t10ω1θ  (tkτ)[2τt1t22Δt2u0+2τt0t2Δt2u1+2τt0t12Δt2u2]dτ   +k1j=1tj+1tjω1θ  (tkτ)[2τtjtj+12Δt2uj1+2τtj1tj+1Δt2uj+2τtj1tj2Δt2uj+1]dτ=1Δt2t10ω1θ  (tkτ)[(τt12)Δu2(τt32)Δu1]dτ   +1Δt2k1j=1tj+1tjω1θ  (tkτ)[(τtj12)Δuj+1(τtj+12)Δuj]dτ=1Δt2[t20ω1θ  (tkτ)(τt32)dτ]Δu1+1Δt2[t20ω1θ  (tkτ)(τt12)dτ]Δu2   +1Δt2k1j=2tj+1tjω1θ  (tkτ)[(τtj12)Δuj+1(τtj+12)Δuj]dτ 0DθΔtuk,

    we can conclude that

    0DθΔtuk=Akk1Δu1+Akk2Δu2+k1j=3AkkjΔuj+Ak0Δuk,3kK, (2.8)

    where

    Ak0=1Δt2tktk1ω1θ  (tkτ)(τtk32)dτ,  Akk1=1Δt2t20ω1θ  (tkτ)(τt32)dτ,Akk2=1Δt2[t20ω1θ  (tkτ)(τt12)dτt3t2ω1θ  (tkτ)(τt52)dτ],Akkj=1Δt2[tjtj1ω1θ  (tkτ)(τtj32)dτtj+1tjω1θ  (tkτ)(τtj+12)dτ],j=3,,k1. (2.9)

    According to approximations (2.5), (2.6) and (2.8), the higher order numerical scheme of Eq (2.1) corresponding to the initial value (2.2) is as follows

    0DθΔtuk=fk,k=1,2,,K. (2.10)

    Remark 2.1. The numerical algorithm in this paper is different from [15]. In [15], in order to obtain (3θ)-order numerical scheme for the linear form of (2.1), they changed (2.1) into the following equation

    1Γ(1θ)t0u(s)(ts)θds=u(0)t1θ  Γ(1θ)(1θ)+t0f(s)ds. (2.11)

    And using the linear Lagrange interpolation in [t0,t1] to u(s), they obtain a (3θ)-order numerical scheme. In this paper, we use the Quadratic Lagrange interpolation in [t0,t2] to u(s) in (2.1) and directly obtain the (3θ)-order numerical scheme (2.5) and (2.6).

    Next, we will introduce a lemma to prove the sign of Akkj and judge its size relation.

    Lemma 2.1. Ak0>0, the sign of Ak1 is uncertain, Ak2>Ak3>>Akk1>0.

    Proof. Through calculation, we can get the following results

    Ak0=1Δt2Γ(2θ)tktk1(τtk32)d(tkτ)1θ  =1ΔtθΓ(3θ)(2θ2)>0, (2.12)
    Ak1=1Δt2[tk1tk2ω1θ  (tkτ)(τtk52)dτtktk1ω1θ  (tkτ)(τtk12)dτ]=1ΔtθΓ(3θ)[(6θ)2θ4+θ]1ΔtθΓ(3θ)M(θ). (2.13)

    The sign of Ak1 is determined by the sign of M(θ), which satisfies

    M(θ)=2θ(log2)[2+(6θ)log2]>0,  M(1)=12(15log2)<0.

    Therefore M(θ) is a decrease function on θ(0,1). By M(0)=2>0 and M(1)=12<0, we know that M(θ) has only one zero θ0 in (0,1), and M(θ)>0 if θ(0,θ0) and M(θ)<0 if θ(θ0,1), which agrees with the conclusion of Lemma 2.1.

    For j=3,,k2, we have

    Akkj=1ΔtθΓ(3θ){12(2θ)[(kj1)1θ  2(kj)1θ  +(kj+1)1θ  ]+(kj1)2θ2(kj)2θ+(kj+1)2θ}, (2.14)

    let ˉj=kj, ˉj2, then using Taylor expansion, we can get the following results

    Akkj=1ΔtθΓ(2θ)ˉjθ+i=0aiˉji, (2.15)

    where ai=Πin=0(1θn)1(i+2)![(1)i(i2)+i+42]. Next, we will prove ˉjθ+i=0aiˉji is decrease function on ˉj. Because +i=0aiˉji is a convergence series, and the radius of convergence exists, it can exchange the derivative and the sum.

    (ˉjθ+i=0aiˉji)=ˉjθ[θˉj1+i=0aiˉji++i=1iaiˉji1]=ˉjθ[θˉj1a0+(1+θ)ˉj2a1++i=2(i+θ)aiˉji1].

    We can find +i=2(i+θ)aiˉji1 is an alternating series with positive first term. So +i=2(i+θ)aiˉji1>0, we have

    θˉj1a0+(1+θ)ˉj2a1++i=2(i+θ)aiˉji1>θˉj1a0+(1+θ)ˉj2a1=θˉj1(1θ)+(1+θ)ˉj2(12)(1θ)θ=(1θ)θ[ˉj112(1+θ)ˉj2]>0.

    Therefore, we have already proved ˉjθ+i=0aiˉji is decrease function on ˉj. According to (2.15), we obtain

    Ak2>Ak3>>Akk3. (2.16)

    Next, let's prove Akk3>Akk2>Akk1>0. Through careful calculation, we can get

    Akk2Akk1=1Γ(3θ)ΔtθF1,Akk3Akk2=1Γ(3θ)ΔtθF2,

    where

    F1=32(2θ)(k2)1θ  +12(2θ)(k3)1θ  2(2θ)k1θ  +(k3)2θ3(k2)2θ+2k2θ,F2=12(2θ)k1θ  +32(2θ)(k2)1θ  32(2θ)(k3)1θ  +12(2θ)(k4)1θ  k2θ+3(k2)2θ3(k3)2θ+(k4)2θ,

    According to (1) in Lemma A.1, we get F1>0,F2>0. Therefore, we have

    Akk2Akk1>0,   Akk3Akk2>0. (2.17)

    The following step, we will prove Akk1>0. By carefully calculation, we have

    Akk1=1ΔtθΓ(3θ){2θ2[(k2)1θ  +3k1θ  ]+(k2)2θk2θ}, (2.18)

    According to (2.33), we have

    Akk1>0. (2.19)

    Combining (2.12), (2.13), (2.16), (2.17) with (2.19), we already proved Lemma 2.1.

    We let Δˉuj=Δ(ujρuj1), Δˉu0=Δu0. Let ρ=12(1Ak1Ak0)=12(3+θ64θ21θ  ). Therefore, from (2.8), we can get that

    0DθΔtuk=Ak0Δuk+Ak1Δuk1+Ak2Δuk2++Akk2Δu2+Akk1Δu1=Ak0Δ(ukρuk1)+(Ak1+ρAk0)Δ(uk1ρuk2)++(Akkj+ρAkkj1++ρkjAk0)Δ(ujρuj1)++(Akk1+ρAkk2++ρk1Ak0)Δ(u1ρu0)+(ρAkk1+ρ2Akk2++ρkAk0)Δu0=kj=1ˉAkkjΔˉuj+Δu0ni=1ρiAkki, (2.20)

    where

    ˉAkkj=ki=jρijAkki. (2.21)

    Next, we discuss the relationship between the size of ˉAkkj.

    Lemma 2.2. ˉAk0>ˉAk1>ˉAk2>>ˉAkk1>0.

    Proof. When k=3, we only need to prove ˉA30>ˉA31>ˉA32>0. By careful calculation, we can get

    ˉA31ˉA30=12ΔtθΓ(3θ)[32(4θ)(4+θ)31θ  +(6θ)2θ].

    According to (2) in Lemma A.1, we get

    ˉA30>ˉA31. (2.22)

    By direct calculation, we have

    ˉA32ˉA31=A32A31+ρ(ˉA31ˉA30)=4θ2ΔtθΓ(3θ)[34+5θ42(4θ)31θ  +(4+θ)(6θ)(4θ)231θ  2θ(6θ)2(4θ)2(2θ)2].

    According to (3) in Lemma A.1, we have

    ˉA31>ˉA32. (2.23)

    Next, we will prove ˉA32=A32+ρˉA31>0. According to (2.41), we get ˉA31>0. Hence, we just need to prove A32>0. By calculation, we have A32=12ΔtθΓ(3θ)[4θθ32θ]>0, therefore,

    ˉA32>0. (2.24)

    Combining (2.22), (2.23) and (2.24), when k=3, we know that Lemma 2.2 holds.

    When k4, through careful calculation, we can draw the following conclusions

    ˉAk1ˉAk0=2θ2ΔtθΓ(3θ)[32+(θ6)2θ4θ]=2θ2ΔtθΓ(3θ)ρ. (2.25)

    According to (4) in Lemma A.1, we can get ρ>0, we have ˉAk0>ˉAk1. Because

    ˉAk2ˉAk1=Ak2Ak1+ρ(ˉAk1ˉAk0)=18(4θ)ΔtθΓ(3θ)˜f(θ),

    where ˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ  4(4θ)(8θ)31θ  +(6θ)241θ  . According to (5) in Lemma A.1, we can get ˜f(θ)>0. It can be concluded that

    ˉAk1>ˉAk2. (2.26)

    By careful calculation, we can get ˉAk3ˉAk2=Ak3Ak2+ρ(ˉAk2ˉAk1). According to ρ>0, (2.26) and Lemma 2.1, we can get

    ˉAk2>ˉAk3. (2.27)

    By mathematical induction, and Lemma 2.1 we can conclude that

    ˉAk1jˉAk2j=Akk1jAkk2j+ρ(ˉAkk2jˉAkk3j)>0,j=0,,k4.

    So we have

    ˉAk0>ˉAk1>ˉAk2>>ˉAkk1. (2.28)

    Next we prove that ˉAkk1>0. According to (2.21), we have

    ˉAkk1=Akk1+k2i=2ρijAkki+ρk2ˉAk1.

    According to (2.37), we get ˉAk1>0. Using Lemma 2.1, ρ>0, we have

    ˉAkk1>0. (2.29)

    Combining (2.28) with (2.29), we already proved Lemma 2.2.

    Lemma 2.3. There is a constant πA6 such that the discrete kernel satisfies the lower bound.

    ˉAkj1πAΔttkjtkj1ω1θ  (tkτ)dτ=1πAΔtθΓ(2θ)[(j+1)1θ  j1θ  ]. (2.30)

    Proof. When k4, According to (2.21), Lemma 2.1 and Lemma 2.2, we can get

    ˉAkkj=Akkj+k2i=j+1Akkiρij+ρk1jˉAk1Akkj.

    For j=1, we have

    ˉAkk1Akk1=1ΔtθΓ(3θ){2θ2[(k2)1θ  +3k1θ  ]+(k2)2θk2θ}, (2.31)

    let ˜j=k1,˜j3, using Taylor formula to expand the calculation, we can get the following results

    Akk1=˜jθΔtθΓ(2θ)+i=0ai, (2.32)

    where ai=Πii=0(1θi)1(i+2)!(1˜j)i[(1)i+1i2+3i+42], a0=(1θ), a1=8(1θ)(θ)12˜j, +i=2ai is a convergent alternating series, and a2>0, therefore, we have 0+i=2aia2, so

    Akk1=˜jθΔtθΓ(2θ)(a0+a1++i=2ai)˜jθΔtθΓ(1θ)(1836)=˜jθΔtθΓ(1θ)97. (2.33)

    We have

    πA97. (2.34)

    For j=2, we have

    ˉAkk2Akk2=1ΔtθΓ(3θ){2θ2[(k3)1θ  2(k2)1θ  k1θ  ]+(k3)2θ2(k2)2θ+k2θ},

    let ˆj=k2,ˆj2, and then using Taylor expansion, similar to j=1, it can be obtained by calculation

    Akk2=ˆjθΔtθΓ(3θ)+i=0ai,

    where ai=Πin=0(1θn)1(i+2)!(1ˆj)i[i2(1)i+2i22i+1], because +i=3ai is a convergent staggered series, and a3>0, 0<+i=3ai<a3, we get

    ˉAkk2ˆjθΔtθΓ(2θ)[a0+a1+a2]ˆjθΔtθΓ(2θ)[151221242]=ˆjθΔtθΓ(2θ)3748,

    therefore, we have

    πA4837. (2.35)

    For j=3,4,,k2, we have

    ˉAkkjAkkj=1ΔtθΓ(3θ){12(2θ)[(kj1)1θ  2(kj)1θ  +(kj+1)1θ  ]+(kj1)2θ2(kj)2θ+(kj+1)2θ},

    let ˉj=kj, ˉj2, we have

    Akkj=ˉj1θ  ΔtθΓ(3θ){12(2θ)[(11ˉj)1θ  2+(1+1ˉj)1θ  ]+ˉj[(11ˉj)2θ2+(1+1ˉj)2θ]},

    then using Taylor expansion, we can get the following results

    Akkj=ˉjθΔtθΓ(2θ)+i=0ai,

    where ai=Πin=0(1θn)1(i+2)!(1ˉj)i[(1)i(i2)+i+42]. Because +i=2ai is a convergent alternating series and a2>0,0+i=2aia2, so

    ˉAkkjˉjθΔtθΓ(1θ)11θ  [1θ+33!(1θ)(θ)1ˉj]ˉjθΔtθΓ(1θ)(114),

    so we can get

    πA43. (2.36)

    When j=k1,

    ˉAk1=ρˉAk0+Ak1=1ΔtθΓ(3θ)[θ41+(3θ2)2θ]12ΔtθΓ(1θ)[θ41+(3θ2)2θ]14ΔtθΓ(1θ). (2.37)

    Therefore, we have

    πA4. (2.38)

    When j=k, according to (2.30), we can get

    ˉAk0=Ak0=1ΔtθΓ(3θ)(2θ2)12ΔtθΓ(2θ)(2θ2)1ΔtθΓ(2θ). (2.39)

    We have

    πA1. (2.40)

    When k=3, ˉA30 is in (2.39), so 1πA1. According to (2.30), we only need

    ˉA31=ρˉA30+A31=1ΔtθΓ(3θ)[θ41+(θ23)2θ+(θ2+2)31θ  ]=21θ  1ΔtθΓ(2θ)θ41+(θ23)2θ+(θ2+2)31θ  (21θ  1)(2θ)21θ  12ΔtθΓ(2θ)1(21θ  1)(2θ)21θ  14ΔtθΓ(2θ). (2.41)

    Therefore, we have

    πA4. (2.42)

    Next, we calculate ˉA32,

    ˉA32=A32+ρˉA31=1ΔtθΓ(3θ)18(4θ)[(θ4)2+4(4θ)(θ6)2θ+4(6θ)2(2θ)2   +[6(4θ)24(6θ)(4+θ)2θ]31θ  ]31θ  21θ  ΔtθΓ(2θ)18(4θ)(31θ  21θ  )(2θ)(2+θ)   [θ36θ2+32+4(θ3+17θ276θ+96)2θ+4(θ34θ272)22θ]31θ  21θ  ΔtθΓ(2θ)116(4θ)(2+θ)[θ36θ2+32+4(θ3+17θ276θ+96)2θ   +4(θ34θ272)22θ]31θ  21θ  ΔtθΓ(2θ)2416(4θ)(2+θ)31θ  21θ  6ΔtθΓ(2θ).

    Therefore, we have

    πA6. (2.43)

    Combining (2.33)–(2.38) with (2.12)–(2.43), we have already proved Lemma 2.3.

    Now we turn to derive an estimate for the truncation errors of the scheme (2.10). We start with deriving an error estimate for the finite difference operator 0DθΔtuk.

    Theorem 3.1. Assume u(t)C3[0,T]. Let

    rk(Δt)=0Dθtu(tk)0DθΔtuk,k1. (3.1)

    Then there exists a constant Cu depending only on the function u, such that for all Δt>0,

    |rk(Δt)|CuΔt3θ. (3.2)

    Proof. Our error estimation will be established on the following Taylor theorem

    u(t)J[t0,t2]u(t)=u(3)(ξ(τ))6(τt0)(τt1)(τt2),τ[t0,t2], (3.3)

    where ξ(τ) is a function defined on [t0,t2] with range (t0,t2). We first estimate r1(Δt)

    |r1(Δt)|=|1Γ(1θ){t10u(τ)(t1τ)θdτt10[J[t0,t2]u(τ)](t1τ)θdτ}|=θΓ(1θ)|t10[u(τ)J[t0,t2]u(τ)](t1τ)θ1dτ|=θΓ(1θ)|t10u(3)(ξ(τ))6(τt0)(τt2)(t1τ)θdτ|θΓ(1θ)maxt[0,T]|u(3)(t)|6t10τ(t2τ)(t1τ)θdτ=θ3(3θ)Γ(2θ)maxt[0,T]|u(3)(t)|Δt3θCuΔt3θ.

    This proves (3.2) for k=1. The case k=2 can be similarly proven, and here we omit the details.

    When k3, we have

    |rk(Δt)|=|0Dθtu(tk)0DθΔtu(tk)|=1Γ(1θ)|t10(tkτ)θ[u(τ)J[t0,t2]u(τ)]dτ+k1j=1tj+1tj(tkτ)θ[u(τ)J[tj1,tj+1]u(τ)]dτ||M+N|.

    Similar to the proof of |r1(Δt)|, we have

    |M|CuΔt3θ. (3.4)

    For N, we have

    |N|=1Γ(1θ)|k1j=1tj+1tj(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]|=1Γ(1θ)|k2j=1tj+1tj(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]+tktk1(tkτ)θd[u(τ)J[tj1,tj+1]u(τ)]|=1Γ(1θ){|k2j=1tj+1tju(3)(ξ)6(τtj1)(τtj)(τtj+1)d(tkτ)θtktk1u(3)(ξ)6(τtk2)(τtk1)(τtk)d(tkτ)θ|}3θ27Γ(1θ)maxt[0,T]|u(3)(t)|Δt3tj+1tj(tkτ)θ1dτ+θ6Γ(1θ)maxt[0,T]|u(3)(t)||Δt2tktk1(tkτ)θdτ3θ27Γ(1θ)maxt[0,T]|u(3)(t)|Δt3θ+θ6Γ(1θ)maxt[0,T]|u(3)(t)|Δt3θ.

    In the above inequality, we use |(τtj1)(τtj)(τtj+1)|239Δt3.

    We can get the following conclusions

    |rk(Δt)|=|M+N|CuΔt3θ. (3.5)

    Complete the proof of Theorem 3.1.

    In order to obtain the convergence analysis, we now introduce an important tool: complementary discrete convolution kernel. Because of ωθωβ=ωθ+β, therefore

    tsωθ(tμ)ω1θ  (μs)dμ=ω1(ts)=1,for0<s<t<. (4.1)

    A class of complementary discrete convolution kernels Pkki with the same properties are given

    ki=mPkkiˉAiim1,for3mkK. (4.2)

    According to [3], we have

    Pk0=1ˉAk0,Pki=1ˉAki0i1j=0(ˉAkjij1ˉAkjij)Pkj,for1ik3, (4.3)

    where Pki>0,i=0,1,...,k3 was obtained from Lemma 2.2.

    Next, we introduce Lemma 2.1 of [3], as follows.

    Lemma 4.1. In (4.3), the discrete kernel Pki has the property of (4.2) and satisfies the following conditions

    0PkkiπAΓ(2θ)Δtθ,1ikN, (4.4)
    kj=1Pkkiω1θ  (ti)πA,1kK. (4.5)

    Now we can analyze the convergence of 0DθΔtuj. First, we consider f(t,u) and assume that it satisfies the following differential mean value theorem and that there exists a constant L, such that

    f(tk,u(tk))f(tk,uk)=Lk(u(tk)uk)=Lkek,k=1,2,...,K, (4.6)

    and |Lk|L,k1.

    Theorem 4.1. Assumes that u is the exact solution of (2.1) and (2.2), and {uk}Kk=1 is the numerical solution of (2.10). If θ(0,1) and step Δt satisfy

    Δt1θ2πAΓ(2θ)Λ, (4.7)

    where πA6, and Λ=12L, then the following error estimates hold,

    |u(tk)uk|CΔt3θ, (4.8)

    here C is only dependent on the function u.

    Proof. Because Δu1=u1u0,Δu2=u2u1, through calculation, we can get

    {0Dθtu(t1)=B11u0+(B11B21)u1+B21u2,0Dθtu(t2)=B12u0+(B12B22)u1+B22u2, (4.9)

    let B11=ΔtθˉD10, B11B21=ΔtθˉD11, B21=ΔtθˉD12; B12=ΔtθD10, B12B22=ΔtθD11, B22=ΔtθD12. It is easy to obtain by calculation

    {ΔtθˉD11e1+ΔtθˉD12e2=f(t1,u(t1))f(t1,u1)r1(Δt),ΔtθD11e1+ΔtθD12e2=f(t2,u(t2))f(t2,u2)r2(Δt). (4.10)

    Among (4.6), then we can get

    {ΔtθˉD11e1+ΔtθˉD12e2=L1e1r1(Δt),ΔtθD11e1+ΔtθD12e2=L2e2r2(Δt).

    That is

    (e1e2)=ΔtθD1(L1e1r1(Δt)L2e2r2(Δt))ΔtθD1L1e1r1(Δt)L2e2r2(Δt)ΔtθD1(Lmax{|e1|,|e2|}+max{|r1(Δt)|,|r2(Δt)|}), (4.11)

    where D=(ˉD11ˉD12D11D12). Since ˉD11D12ˉD12D11 is bounded, then

    max{|e1|,|e2|}C1Δtθmax{|r1(Δt)|,|r2(Δt)|}, (4.12)

    therefore,

    {|e1|CΔt3θ,|e2|CΔt3θ, (4.13)

    when j=1,2, (4.8) is proved.

    When i3, let's set ˉei=eiρei1,i=1,...,j, and ˉe0=e0 for (4.13), we have

    {|ˉe1|=|e1ρe0|CΔt3θ,|ˉe2|=|e2ρe1|53CΔt3θ. (4.14)

    Combining with (2.20) and (4.6), ˉej, j3, satisfy

    kj=1ˉAkkjΔˉej=f(tk,u(tk))f(tk,uk)rk(Δt). (4.15)

    Because of ej=jn=1ρjnˉen, then

    kj=3ˉAkkjΔˉej=Lkkj=3ρkjˉejˉrk(Δt), (4.16)

    where

    ˉrk(Δt)=rk(Δt)Lk2j=1ρkjˉej+2j=1ˉAkkjΔˉej, (4.17)

    by multiplying 2ˉek on both sides of (4.16), from Lemma A.1 in Liao [3], we can obtain that

    2ˉekkj=3ˉAkkjΔˉejkj=3ˉAkkjΔ(|ˉej|2),2ˉek[Lkkj=3ρkjˉejˉrk(Δt)]Lkj=3ρkj(|ˉej|2+|ˉek|2)+2|ˉek||ˉrk(Δt)|4Lkj=3ρkj|ˉej|2+2|ˉek||ˉrk(Δt)|.

    To sum up, we can get the following conclusions

    kj=3ˉAkkjΔ(|ˉej|2)4Lkj=3ρkj|ˉej|2+2|ˉek||ˉrk(Δt)|. (4.18)

    We change k to i, and then multiply both sides of (4.18) by a ki=3Pkki, then

    ki=3Pkkiij=3ˉAiijΔ(|ˉej|2)ki=3Pkkiij=34Lρij|ˉej|2+2ki=3Pkki|ˉei||ˉri(Δt)|. (4.19)

    From (4.2) and Δuj=ujuj1, we can get the following results

    ki=3Pkkiij=3ˉAiijΔ(|ˉej|2)=kj=3Δ(|ˉej|2)ki=jPkkiˉAiij=kj=3Δ(|ˉej|2)=|ˉek|2|ˉe2|2, (4.20)

    substituting (4.20) into (4.19) yields

    |ˉek|24Lki=3Pkkiij=3ρij|ˉej|2+|e2|2+2ki=3Pkki|ˉei||ˉri(Δt)|. (4.21)

    Next, we will prove it by mathematical induction

    |ˉek|2Eθ(2πAΛtθk)(|ˉe2|+2max3jkki=3pjji|ˉri(Δt)|)FkGk,fork2, (4.22)

    where Λ=12L, Eθ stands for the Mittag-Leffler function, we usually define it as: Eθ(z)=j=0ZjΓ(1+jθ), and Eθ(0)=1,ZR, and Z>0,Eθ(Z)>0, so FkFk12 for all k2. where Gk=|ˉe2|+2max3jkki=3Pjji|ˉri(Δt)|.

    When k=3, (i) if |ˉe3||e2|,G3=|e2|+2P30|ˉr3(Δt)||ˉe2|, we have |ˉe3||ˉe2|G3F3G3,

    (ii) if |ˉe3|>|e2|, from (4.21), it can be concluded that

    |ˉe3|24LP30|ˉe3|2+|ˉe2|2+2P30|ˉe3||ˉr3(Δt)|.

    According to (4.4), it can be concluded that

    4LP30πAΓ(2θ)ΔtθΛ12Δt1θ2πAΓ(2θ)Λ,

    so |ˉe3|2(|ˉe2|+2P30|ˉr3(Δt)|)2G3F3G3.

    When 4kK, we assume that: |ˉej|FjGj, for 3jk1, let |ˉej(k)|=max2ik1|ˉei|.

    (i) When |ˉek||ˉej(k)|, because Fj and Gj monotonically increase with j, so

    |ˉek||ˉej(k)|Fj(k)Gj(k)FkGk.

    (ii) When |ˉek||ˉej(k)|, from (4.21) we can obtain that

    |ˉek|2|ˉek|k1i=3Pkkiij=34Lρij|ˉej|+|ˉek|2Pk0kj=34Lρij+|ˉek||ˉe2|+2|ˉek|ki=3Pkki|ˉri(Δt)|.

    From (4.4), we give the limit of the maximum step size, which means that

    Pk0kj=34LρijπAΓ(2θ)ΔtθΛ<12Δt1θ2πAΓ(2θ)Λ.

    Using of Lemma 2.3 in [3] for any real μ>0,

    k1i=3PkkiEθ(μtθi)πAEθ(μtθk)1μ,for3kK.

    We have

    |ˉek|2k1i=3Pkkiij=34Lρij|ˉej|+2|ˉe2|+4ki=3Pkki|ˉri(Δt)|2Λk1i=3PkkiFiGk+2Gk2k1i=3Pkkiij=34LθijFjGj+2Gk[4Λk1i=3PkkiEθ(2πAΛtθi)+2]Gk{4ΛπA2πAΛ[Eθ(2πAtθk)1]+2}Gk=2Eθ(2πAΛtθk)Gk=FkGk.

    To sum up, the proof of (4.22) is complete.

    Next, we carefully estimate |ˉri(Δt)|. According to (4.17), we have

    |ˉrk(Δt)||rk(Δt)|+L2j=1ρkj|ˉej|+(ˉAkk1+ˉAkk2)|ˉe1|+ˉAkk2|ˉe2|. (4.23)

    Next, we estimate the upper bounds of ˉA(k)k1 and ˉA(k)k2. According to (2.37), Lemma 2.1 and (4) in Lemma A.1, we can obtain

    ˉAkk2=k2j=2Akkjρj2+ρk3ˉAk1Ak2k2j=2ρj2+ρk3ˉAk13Ak2+1ΔtθΓ(3θ)[θ41+(3θ2)2θ]=1ΔtθΓ(3θ)[54(4θ)+92(8θ)3θ+112(6+θ)2θ]1ΔtθΓ(3θ)[5+928+112(6+1)12]=1ΔtθΓ(3θ)1094.

    Therefore, ˉAkk21094ΔtθΓ(3θ). According to Lemma 2.2, we can get ˉAkk1<ˉAkk21094ΔtθΓ(3θ). Therefore, we can obtain from (4.23),

    |ˉri(Δt)|CΔt3θ+5LCΔt3+1092Γ(3θ)Δt3θ+54512Γ(3θ)CΔt3θ[C+5LCΔtθ+100Γ(3θ)C]Δt3θ.

    In the expression on the right side of (4.22),

    ki=3Pjji|ˉri(Δt)|=ki=3Pjjiω1θ  (ti)max3ik|ˉri(Δt)|ω1θ  (ti), (4.24)

    according to (4.5) in Lemma 4.1, we can obtain that

    ki=3Pjjiω1θ  (ti)ki=1Pjjiω1θ  (ti)πA,1kK.

    Because of 1ω1θ  (ti)=Γ(1θ)tθi, therefore (4.24) becomes

    ki=3Pjji|ˉri(Δt)|πAmax3ik|ˉri(Δt)ω1θ  (ti)πA[C+5LCΔtθ+100Γ(3θ)C]Δt3θmax3ikΓ(1θ)tθi. (4.25)

    Based on (4.13), (4.25) and (4.22), we obtain

    |ˉek|2Eθ(2πAΛtθk){53CΔtθ+2πA[C+5LCΔtθ+100Γ(3θ)C]Γ(1θ)tθk}Δt3θ.

    According to ek=ˉekρek1=kn=0ρknˉen, we have |ek|3max0nk|ˉen|. Then the proof is completed.

    In this section, two examples are presented to verify the effectiveness of our numerical method (2.10).

    Example 5.1. we consider the problem (2.1) with u(0)=0, and

    f(t,u(t))=Γ(4+θ)6t3+t3+θu(t), (5.1)
    f(t,u(t))=Γ(4+θ)6t3+t6+2θu2(t). (5.2)

    The exact solution of Eq (2.1) is u(t)=t3+θ. We take T=1 and choose the step size to be Δt=2l,l=7,8,,11. The error we will display is defined by eΔt=maxk=1,2,,K|u(tk)uk|,K=T/Δt.

    In (5.1), f is a linear case of u. From Table 1, the convergence order is computed by log2(e2Δt/eΔt). it is observed that for θ=0.3,0.6 and 0.9, the convergence rates are close to 2.7,2.4 and 2.1, respectively. This is in a good agreement with the theoretical prediction. In Table 2, we can take θ=0.01 and 0.99 and obtain that the convergence rates are close to 2.99 and 2.01. That tells us that as θ0 or 1, the convergence rate still 3θ. In (5.2), f is the nonlinear case of u. From Tables 3 and 4, once again these results confirm that the convergence of the numerical solution is close to of order 3θ for 0<θ<1.

    Table 1.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.1).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 5.1545E-6 - 5.3113E-5 - 3.8546E-4 -
    1256 8.1757E-7 2.6564 1.0192E-5 2.3815 9.0841E-5 2.0851
    1512 1.2874E-7 2.6668 1.9457E-6 2.3891 2.1297E-5 2.0926
    11024 2.0164E-8 2.6746 3.7032E-7 2.3934 4.9804E-6 2.0963
    12048 3.1714E-9 2.6686 7.0366E-8 2.3958 1.1632E-6 2.0981

     | Show Table
    DownLoad: CSV
    Table 2.  Maximum errors eΔt and decay rate as Δt=0.01,0.99 with the right hand side (5.1).
    Δt θ=0.01 Rate θ=0.99 Rate
    1128 3.7469E-8 - 6.6342E-4 -
    1256 5.1473E-9 2.8638 1.6645E-4 1.9947
    1512 7.0275E-10 2.8727 4.1540E-5 2.0025
    11024 8.9656E-11 2.9705 1.0339E-5 2.0063
    12048 1.1231E-11 2.9968 2.5703E-6 2.0081

     | Show Table
    DownLoad: CSV
    Table 3.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.2).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 1.0154E-6 - 1.0167E-5 - 9.1822E-5 -
    1256 1.4941E-7 2.7647 1.9494E-6 2.3828 2.1549E-5 2.0911
    1512 2.3628E-8 2.6607 3.7176E-7 2.3906 5.0416E-6 2.0957
    11024 3.7049E-9 2.6729 7.0707E-8 2.3944 1.1777E-6 2.0978
    12048 5.6279E-10 2.7187 1.3497E-8 2.3891 2.7492E-7 2.0989

     | Show Table
    DownLoad: CSV
    Table 4.  Maximum errors eΔt and decay rate as Δt=0.01,0.99 with the right hand side (5.2).
    Δt θ=0.01 Rate θ=0.99 Rate
    1128 3.2116E-6 - 1.7072E-4 -
    1256 3.9867E-7 3.0099 4.2645E-5 2.0012
    1512 4.9490E-8 3.0099 1.0618E-5 2.0058
    11024 6.1435E-9 3.0099 2.6400E-6 2.0079
    12048 7.6264E-10 3.0099 6.5589E-7 2.0090

     | Show Table
    DownLoad: CSV

    Example 5.2. we consider the problem (2.1) and (2.2) with the following right hand side function

    f(t,u(t))=tθSin1,1α(t)+sin(t)u(t), (5.3)
    f(t,u(t))=tθSin1,1α(t)+sin2(t)u2(t). (5.4)

    The corresponding exact solution is u(t)=sin(t), and Sinα,β(t)=k=1(1)k+1t2k1Γ(α(2k1)+β). f is a linear and nonlinear case of u in (5.3) and (5.4), respectively.

    We repeat the same calculation as in Example 5.1 using the proposed numerical scheme. Tables 5 and 6 show the maximum errors and decay rates of the step size for several a ranging from 0.3 to 0.9. This is consistent with the theoretical prediction. The convergence rate is close to 3θ for 0<θ<1.

    Table 5.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.3).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 4.4720E-7 - 3.3763E-6 - 2.1961E-5 -
    1256 7.1175E-8 2.6515 6.4755E-7 2.3823 5.1721E-6 2.0861
    1512 1.1232E-8 2.6636 1.2353E-7 2.3901 1.2121E-6 2.0931
    11024 1.7643E-9 2.6704 2.3497E-8 2.3942 2.8338E-7 2.0967
    12048 2.6806E-10 2.7185 4.4390E-9 2.4041 6.6105E-8 2.0999

     | Show Table
    DownLoad: CSV
    Table 6.  Maximum errors eΔt and decay rate as functions of Δt with the right hand side (5.4).
    Δt θ=0.3 Rate θ=0.6 Rate θ=0.9 Rate
    1128 5.2408E-7 - 3.5338E-6 - 2.0939E-5 -
    1256 8.5515E-8 2.6155 6.8794E-7 2.3608 4.9822E-6 2.0713
    1512 1.3708E-8 2.6411 1.3230E-7 2.3784 1.1738E-6 2.0855
    11024 2.1730E-9 2.6573 2.5279E-8 2.3878 2.7517E-7 2.0928
    12048 3.4105E-10 2.6716 4.8110E-9 2.3935 6.4316E-8 2.0970

     | Show Table
    DownLoad: CSV

    In this paper, we presented a high order numerical method for solving Caputo nonlinear fractional ordinary differential equations. The numerical method was constructed by using the Quadratic Lagrange interpolation. By careful error estimation, the proposed scheme is of order 3θ for 0<θ<1. Finally, numerical experiments are carried out to verify the theoretical prediction. In the future, we plan to apply this kind of methods to 3D fractional partial differential equations with time derivatives based the idea of [21] and [22].

    This research was supported by National Natural Science Foundation of China (Grant No. 11901135, 11961009), Foundation of Guizhou Science and Technology Department (Grant No. [2020]1Y015, [2017]1086), Foundation for graduate students of Guizhou Provincial Department of Education(Grant No.YJSCXJH[2020]136). The authors thank the anonymous referees for their valuable suggestions to improve the quality of this work significantly.

    The authors declare that there is no conflict of interests regarding the publication of this article.

    Lemma A.1.

    (1)F1=32(2θ)(k2)1θ+12(2θ)(k3)1θ2(2θ)k1θ+(k3)2θ3(k2)2θ+2k2θ>0,F2=12(2θ)k1θ+32(2θ)(k2)1θ32(2θ)(k3)1θ+12(2θ)(k4)1θk2θ+3(k2)2θ3(k3)2θ+(k4)2θ>0,(2)32(4θ)(4+θ)31θ+(6θ)2θ>0,(3)34+5θ42(4θ)31θ+(4+θ)(6θ)(4θ)231θ2θ(6θ)2(4θ)2(2θ)2>0,(4)0<ρ=12(3+θ64θ21θ)<23,(5)˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ4(4θ)(8θ)31θ+(6θ)241θ>0,
    (1)F1=32(2θ)(k2)1θ+12(2θ)(k3)1θ2(2θ)k1θ+(k3)2θ3(k2)2θ+2k2θ>0,F2=12(2θ)k1θ+32(2θ)(k2)1θ32(2θ)(k3)1θ+12(2θ)(k4)1θk2θ+3(k2)2θ3(k3)2θ+(k4)2θ>0,(2)32(4θ)(4+θ)31θ+(6θ)2θ>0,(3)34+5θ42(4θ)31θ+(4+θ)(6θ)(4θ)231θ2θ(6θ)2(4θ)2(2θ)2>0,(4)0<ρ=12(3+θ64θ21θ)<23,(5)˜f(θ)=3(4θ)2+6(4θ)(6θ)21θ4(4θ)(8θ)31θ+(6θ)241θ>0,

    Proof. (1) Firstly, we will prove F1>0. let k2=x, Taylor formula is used to obtain the following results

    F1=(2θ)(1θ)xθ1+n=0ni=0(θi)(1x)nan,

    where an=12(n+1)(1)n8(n+1)2n(n+3)!=12(n+1)[(1)n162n](n+3)!<0, so +n=0ni=0(θi)(1x)nan is an alternating series with positive first term, so satisfies +n=0ni=0(θi)(1x)nan>0, so F1>0.

    Next, we will prove F2>0. Let k2=ˉx, Taylor formula is used to obtain the following results

    F2=(2θ)(1θ)ˉxθ1+n=0ni=0(θi)(1ˉx)n1(n+3)!bn,

    where bn=32(n+1)(1)n+1+2(n1)2n[1+(1)n],b0=112,b1=3, when n2, n is odd, bn>0; n is even, bn=(n1)[42n32]3>[832]3>0, so +n=1ni=0(θi)(1ˉx)n1(n+3)!bn is an alternating series with positive first term. So +n=0ni=0(θi)(1ˉx)n1(n+3)!bn>(θ)13!b0+0>0. We get F2>0.

    (2) Let f(θ)=32(4θ)(4+θ)31θ+(6θ)2θ, through careful calculation, we get

    f(θ)=3231θ[1+(4+θ)(ln3)]+2θ[1+(6θ)(ln2)],

    f(θ)=3θg(θ), where g(θ)=3(ln3)[2(4+θ)ln3]+(32)θ(ln2)[2+(6θ)ln2], because

    g(θ)=3(ln3)2+(32)θ(ln2){[2+(6θ)(ln2)]ln(32)ln2},g(θ)=(32)θ(ln2)ln(32){[2+(6θ)(ln2)]ln(32)2ln2}{32}θ(ln2)ln(32)H(θ),

    We can get H(θ)=ln2ln(32)<0, so we get H(θ)>H(1)>0, therefore g(θ)<g(1)<0, g(θ) monotone decreasing, g(1)<g(θ)<g(0)=3.6227<0, that is f(θ) monotone decreasing, 0<f(1)<f(θ)<f(0), where f(1)=3+5ln352ln2>0, that is f(θ) monotone increasing, f(0)<f(θ)<f(1), where f(0)=0, to sum up, we can get f(θ)>0.

    (3) Let f1(θ)=34+5θ42(4θ)31θ+(4+θ)(6θ)(4θ)231θ2θ(6θ)2(4θ)2(2θ)2-34+a131θ+a231θ2θ+a3(2θ)2, where

    a1=5θ42(4θ),a2=(4+θ)(6θ)(4θ)2,a3=(6θ)2(4θ)2.

    Next, we just need to prove that f1(θ)>0, using a Taylor expansion yields

    f1(θ)=34+[a121θ+a2212θ](1+12)1θ+a3(2θ)2=34+a121θ+a2212θ+a322θ+[a121θ+a2212θ][1θ1!12+(1θ)(θ)2!(12)2+(1θ)(θ)(θ1)3!(12)3+],

    where

    a121θ+a2212θ=2θ1(4θ)2[(5θ4)(4θ)+2(4+θ)(6θ)2θ]2θ1(4θ)2[(5θ4)(4θ)+(4+θ)(6θ)]0.

    Because of 1θ1!12+(1θ)(θ)2!(12)2+(1θ)(θ)(θ1)3!(12)3++k=0ak is an alternating series with positive first term, and +k=0ak=a0+a1++k=2ak, where +k=2ak is an alternating series with positive first term, so 0<+k=2ak<a2, we have

    f1(θ)>34+a322θ+[a121θ+a2212θ][1+1θ1!12+(1θ)(θ)2!(12)2]=34+2θ8(4θ)2[5θ4+49θ3196θ2+368θ192+(θ4+7θ32θ248θ+144)21θ]34+2θ8(4θ)2[f2(θ)+f3(θ)].

    We wish prove f1(θ)>0, only need 34+2θ8(4θ)2[f2(θ)+f3(θ)]>0f2(θ)+f3(θ)>348(4θ)22θ=6(4θ)22θ, that is f2(θ)+f3(θ)>6(4θ)22θ6f4(θ), so to prove f1(θ)>0, just prove f2(θ)+f3(θ)6f4(θ)>0, let's remember ˉf(θ)=f2(θ)+f3(θ)6f4(θ), since ˉf(θ) first increases and then decreases, and the two endpoints are ˉf(0)=0,ˉf(1)=16>0, ˉf(θ)>0 is always true, so f1(θ)>0.

    (4) Firstly, we will prove ρ>0. Because ρ=3(4θ)+(θ6)21θ2(4θ), let ˉg(θ)=3(4θ)+(θ6)21θ, by calculation, we can get ˉg(θ) is monotonically increasing function on θ, that is ˉg(θ)>ˉg(0)=0. therefore, we have ρ>0.

    In addition, ρ23=5(4θ)+3(θ6)21θ6(4θ)ˉg1(θ)6(4θ). we can directly find ˉg1(θ)<ˉg1(1)=0, so 0<ρ<23.

    (5) We use Lagrange's mean value theorem 41θ>43+θ31θ, and we have

    ˜f(θ)>3(4θ)2+6(4θ)(6θ)21θ+[4(4θ)(8θ)+4(6θ)23+θ]31θ. (A.1)

    Let a1=3(4θ)2,a2=6(4θ)(6θ),a3=43+θ[(4θ)(8θ)(3+θ)+(6θ)2],a4=1+1θ1!12+(1θ)(θ)2(12)2. (A.1) is becomes

    ˜f(θ)>a1+a221θ+a331θ=a1+21θ[a2+a3(1+12)1θ]=a1+21θ{a2+a3[1+1θ1!12+(1θ)(θ)2(12)2]+a3[(1θ)(θ)(θ1)3!(12)3+(1θ)(θ)(θ1)(θ2)4!(12)4+]}=a1+21θ{a2+a3a4+a3(1θ)θ+k=0Πki=0(θ1i)1(k+3)!(12)k+3}a1+21θ{a2+a3a4+a3(1θ)θ+k=0bk},

    where bk=Πki=0(θ1i)1(k+3)!(12)k+3, and b0=(θ1)3!(12)3>0, |bk+1bk|=θ+2+k2(k+4)<1, so +k=0bk is a convergent alternating series, that is 0<+k=0bk<b0. Because of a3<0, so we get

    ˜f(θ)>a1+21θ{a2+a3a4+a3(1θ)θb0}=a1+21θ{a2+a3[a4+(1θ)θ(θ1)3!(12)3]}a1+21θ{a2+a3a5}˜f1(θ),

    where a5=a4+(1θ)θ(θ1)3!(12)3=48+24(1θ)6(θθ2)+(θθ3)48, by careful calculation, we can get

    ˜f1(θ)=3612(3+θ)[θ35θ28θ+48]+21θ112(3+θ)[θ616θ5+97θ4278θ3+88θ2+732θ+864]=112(3+θ){36(θ35θ28θ+48)+21θ(θ616θ5+97θ4278θ3+88θ2+732θ+864)}112(3+θ)˜f2(θ),

    because θ(0,1), so θ3<θ2,θ6>0.

    ˜f2(θ)>36(θ25θ28θ+48)+21θ(016θ5+97θ4278θ3+88θ2+732θ+864)=144(θ2+2θ12)+21θ(16θ5+97θ4278θ3+88θ2+732θ+864)=144(θ2+2θ12)+21θ˜f4(θ)˜f3(θ),

    where ˜f4(θ)=16θ5+97θ4278θ3+88θ2+732θ+864, by careful calculation, we can get ˜f3(θ)<0, that is ˜f3(θ) monotone decreasing, so ˜f3(1)<˜f3(θ)<˜f3(0), by direct calculation, we can get ˜f3(θ)=144(2θ+2)+21θ[˜f4(θ)(ln2)+˜f4(θ)],˜f3(0)>0,˜f3(0)<0, so ˜f3(0) changes from positive to negative, ˜f3(0) increases first and then decreases, ˜f3(0)=0,˜f3(1)=223>0, that is ˜f3(θ)>0 established, so ˜f(θ)>0.



    [1] WHO, Global Tuberculosis Report 2021, 2021. Available from: https://www.who.int/publications/i/item/9789240037021
    [2] WHO, Coronavirus Disease (COVID-19) Pandemic, 2019. Available from: https://www.who.int/emergencies/diseases/novel-coronavirus-2019
    [3] D. Tian, Y. Sun, J. Zhou, Q. Ye, The global epidemic of the SARS-CoV-2 Delta variant, key spike mutations and immune escape, Front. Immunol., 12 (2021), 5001. https://doi.org/10.3389/fimmu.2021.751778 doi: 10.3389/fimmu.2021.751778
    [4] L. Bian, Q. Gao, F. Gao, Q. Wang, Q. He, X. Wu, et al., Impact of the Delta variant on vaccine efficacy and response strategies, Expert Rev. Vaccines, 20 (2021), 1201–1209. https://10.1080/14760584.2021.1976153 doi: 10.1080/14760584.2021.1976153
    [5] S. Shiehzadegan, N. Alaghemand, M. Fox, V. Venketaraman, Analysis of the Delta Variant B.1.617.2 COVID-19, Clin. Pract., 11 (2021), 778–784. https://doi.org/10.3390/clinpract11040093 doi: 10.3390/clinpract11040093
    [6] WHO, Tracking SARS-CoV-2 Variants, 2022. Available from: https://www.who.int/activities/tracking-SARS-CoV-2-variants
    [7] V. T. Hoang, T. L. Dao, P. Gautret, Recurrence of positive SARS-CoV-2 in patients recovered from COVID-19, J. Med. Virol., 92 (2020), 2366–2367.
    [8] F. Zhou, T. Yu, R. Du, G. Fan, Y. Liu, Z. Liu, et al., Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study, Lancet, 395 (2020), 1054–1062. https://10.1016/s0140-6736(20)30566-3
    [9] P. Y. Chia, S. W. X. Ong, C. J. Chiew, L. W. Ang, J. M. Chavatte, T. M. Mak, et al., Virological and serological kinetics of SARS-CoV-2 Delta variant vaccine breakthrough infections: A multicentre cohort study, Clin. Microbiol. Infec., 28 (2021). https://doi.org/10.1016/j.cmi.2021.11.010 doi: 10.1016/j.cmi.2021.11.010
    [10] H. Kato, K. Miyakawa, N. Ohtake, H. Go, Y. Yamaoka, S. Yajima, et al., Antibody titers against the Alpha, Beta, Gamma, and Delta variants of SARS-CoV-2 induced by BNT162b2 vaccination measured using automated chemiluminescent enzyme immunoassay, J. Infect. Chemother., 28 (2022), 273–278. https://doi.org/10.1016/j.jiac.2021.11.021 doi: 10.1016/j.jiac.2021.11.021
    [11] D. S. Khoury, D. Cromer, A. Reynaldi, T. E. Schlub, A. K. Wheatley, J. A. Juno, et al., Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection, Nat. Med., 27 (2021), 1205–1211. https://doi.org/10.1038/s41591-021-01377-8 doi: 10.1038/s41591-021-01377-8
    [12] Z. Fang, Y. Zhang, C. Hang, J. Ai, S. Li, W. Zhang, Comparisons of viral shedding time of SARS-CoV-2 of different samples in ICU and non-ICU patients, J. Infection, 81 (2020), 147–178. https://doi.org/10.1016/j.jinf.2020.03.013 doi: 10.1016/j.jinf.2020.03.013
    [13] Y. Wang, L. Zhang, L. Sang, F. Ye, S. Ruan, B. Zhong, et al., Kinetics of viral load and antibody response in relation to COVID-19 severity, J. Clin. Invest., 130 (2020), 5235–5244. https://doi.org/10.1172/JCI138759 doi: 10.1172/JCI138759
    [14] Y. Liu, L. M. Yan, L. Wan, T. X. Xiang, A. Le, J. M. Liu, et al., Viral dynamics in mild and severe cases of COVID-19, Lancet Infect. Dis., 20 (2020), 656–657. https://doi.org/10.1016/S1473-3099(20)30232-2 doi: 10.1016/S1473-3099(20)30232-2
    [15] Y. Fu, P. Han, R. Zhu, T. Bai, J. Yi, X. Zhao, et al., Risk factors for viral RNA shedding in COVID-19 patients, Eur. Respiratory J., 56 (2020), 2001190. https://doi.org/10.1183/13993003.01190-2020 doi: 10.1183/13993003.01190-2020
    [16] K. Xu, Y. Chen, J. Yuan, P. Yi, C. Ding, W. Wu, et al., Factors Associated With Prolonged Viral RNA Shedding in Patients with Coronavirus Disease 2019 (COVID-19), Clin. Infect. Dis., 71 (2020), 799–806. https://doi.org/10.1093/cid/ciaa351 doi: 10.1093/cid/ciaa351
    [17] A. Mondi, P. Lorenzini, C. Castilletti, R. Gagliardini, E. Lalle, A. Corpolongo, et al., Risk and predictive factors of prolonged viral RNA shedding in upper respiratory specimens in a large cohort of COVID-19 patients admitted to an Italian reference hospital, Int. J. Infect. Dis., 105 (2021), 532–539. https://doi.org/10.1016/j.ijid.2021.02.117 doi: 10.1016/j.ijid.2021.02.117
    [18] H. Long, J. Zhao, H. L. Zeng, Q. B. Lu, L. Q. Fang, Q. Wang, et al., Prolonged viral shedding of SARS-CoV-2 and related factors in symptomatic COVID-19 patients: a prospective study, BMC Infect. Dis., 21 (2021), 1282. https://doi.org/10.1186/s12879-021-07002-w doi: 10.1186/s12879-021-07002-w
    [19] J. Gong, H. Dong, D. K. Wang, F. E. Lu, Z. Y. Huang, K. Fang, et al., Characteristics of Viral Shedding in Respiratory Samples and Specific Antibodies Production in 564 COVID-19 Patients, Curr. Med. Sci., 41 (2021), 46–50. https://doi.org/10.1007/s11596-021-2316-3 doi: 10.1007/s11596-021-2316-3
    [20] R. Ke, P. P. Martinez, R. L. Smith, L. L. Gibson, C. J. Achenbach, S. McFall, et al., Longitudinal analysis of SARS-CoV-2 vaccine breakthrough infections reveal limited infectious virus shedding and restricted tissue distribution, Open Forum Infect. Dis., 9 (2022). https://doi.org/10.1093/ofid/ofac192 doi: 10.1093/ofid/ofac192
    [21] K. Li, S. Luo, Dynamic prediction of Alzheimer's disease progression using features of multiple longitudinal outcomes and time-to-event data, Stat. Med., 38 (2019), 4804–4818. https://doi.org/10.1002/sim.8334 doi: 10.1002/sim.8334
    [22] A. Thomas, G. K. Vishwakarma, A. Bhattacharjee, Joint modeling of longitudinal and time-to-event data on multivariate protein biomarkers, J. Comput. Appl. Math., 381 (2021), 113016. https://doi.org/10.1016/j.cam.2020.113016 doi: 10.1016/j.cam.2020.113016
    [23] D. Rizopoulos, The R Package JMbayes for Fitting Joint Models for Longitudinal and Time-to-Event Data using MCMC, preprient, arXiv: 1404.7625.
    [24] Y. L. Lee, C. H. Liao, P. Y. Liu, C. Y. Cheng, M. Y. Chung, C. E. Liu, et al., Dynamics of anti-SARS-Cov-2 IgM and IgG antibodies among COVID-19 patients, J. Infect., 81 (2020), e55–e58. https://doi.org/10.1016/j.jinf.2020.04.019 doi: 10.1016/j.jinf.2020.04.019
    [25] J. Shang, Y. Wan, C. Luo, G. Ye, Q. Geng, A. Auerbach, et al., Cell entry mechanisms of SARS-CoV-2, Proc. Natl. Acad. Sci. USA, 117 (2020), 11727–11734. https://doi.org/10.1073/pnas.2003138117 doi: 10.1073/pnas.2003138117
    [26] Z. Wang, F. Schmidt, Y. Weisblum, F. Muecksch, C. O. Barnes, S. Finkin, et al., mRNA vaccine-elicited antibodies to SARS-CoV-2 and circulating variants, Nature, 592 (2021), 616–622. https://doi.org/10.1038/s41586-021-03324-6 doi: 10.1038/s41586-021-03324-6
    [27] Y. Li, G. Wang, N. Li, Y. Wang, Q. Zhu, H. Chu, et al., Structural insights into immunoglobulin M, Science, 367 (2020), 1014–1017. https://doi.org/10.1126/science.aaz5425 doi: 10.1126/science.aaz5425
    [28] C. Gaebler, Z. Wang, J. C. C. Lorenzi, F. Muecksch, S. Finkin, M. Tokuyama, et al., Evolution of antibody immunity to SARS-CoV-2, Nature, 591 (2021), 639–644. https://doi.org/10.1038/s41586-021-03207-w doi: 10.1038/s41586-021-03207-w
    [29] J. Carrillo, N. Izquierdo-Useros, C. Ávila-Nieto, E. Pradenas, B. Clotet, J. Blanco, Humoral immune responses and neutralizing antibodies against SARS-CoV-2; implications in pathogenesis and protective immunity, Biochem. Biophs. Res. Commun., 538 (2021), 187–191. https://doi.org/10.1016/j.bbrc.2020.10.108 doi: 10.1016/j.bbrc.2020.10.108
    [30] M. A. Tortorici, M. Beltramello, F. A. Lempp, D. Pinto, H. V. Dang, L. E. Rosen, et al., Ultrapotent human antibodies protect against SARS-CoV-2 challenge via multiple mechanisms, Science, 370 (2020), 950–957. https://doi.org/10.1126/science.abe3354 doi: 10.1126/science.abe3354
    [31] Z. Ku, X. Xie, P. R. Hinton, X. Liu, X. Ye, A. E. Muruato, et al., Nasal delivery of an IgM offers broad protection from SARS-CoV-2 variants, Nature, 595 (2021), 718–723. https://doi.org/10.1038/s41586-021-03673-2 doi: 10.1038/s41586-021-03673-2
    [32] T. Fiolet, Y. Kherabi, C. J. MacDonald, J. Ghosn, N. Peiffer-Smadja, Comparing COVID-19 vaccines for their characteristics, efficacy and effectiveness against SARS-CoV-2 and variants of concern: A narrative review, Clin. Microbiol. Infect., 28 (2022), 202–221. https://doi.org/10.1016/j.cmi.2021.10.005 doi: 10.1016/j.cmi.2021.10.005
    [33] M. Levine-Tiefenbrun, I. Yelin, R. Katz, E. Herzel, Z. Golan, L. Schreiber, et al., Initial report of decreased SARS-CoV-2 viral load after inoculation with the BNT162b2 vaccine, Nat. Med., 27 (2021), 790–792. https://doi.org/10.1038/s41591-021-01316-7 doi: 10.1038/s41591-021-01316-7
  • This article has been cited by:

    1. Xiaopeng Yi, Chongyang Liu, Huey Tyng Cheong, Kok Lay Teo, Song Wang, A third-order numerical method for solving fractional ordinary differential equations, 2024, 9, 2473-6988, 21125, 10.3934/math.20241026
  • Reader Comments
  • © 2023 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(2190) PDF downloads(96) Cited by(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog