Processing math: 83%
Research article Special Issues

Global stability of an HIV infection model with saturated CTL immune response and intracellular delay

  • In this paper, we consider an HIV infection model with saturated infection rate, intracellular delay and saturated cytotoxic T lymphocyte (CTL) immune response. By calculation, we obtain immunity-inactivated reproduction number R0 and immunity-activated reproduction number R1. By analyzing the distribution of roots of the corresponding characteristic equations, we study the local stability of an infection-free equilibrium, an immunity-inactivated equilibrium and an immunity-activated equilibrium of the model. By constructing suitable Lyapunov functionals and using LaSalle's invariance principle, we show that if R0<1, the infection-free equilibrium is globally asymptotically stable; If R1<1<R0, the immunity-inactivated equilibrium is globally asymptotically stable; If R1>1, the immunity-activated equilibrium is globally asymptotically stable. Sensitivity analyses are carried out to show the effects of parameters on the immunity-activated reproduction number R1 and the viral load.

    Citation: Jian Ren, Rui Xu, Liangchen Li. Global stability of an HIV infection model with saturated CTL immune response and intracellular delay[J]. Mathematical Biosciences and Engineering, 2021, 18(1): 57-68. doi: 10.3934/mbe.2021003

    Related Papers:

    [1] Zhiyan Ding, Hichem Hajaiej . On a fractional Schrödinger equation in the presence of harmonic potential. Electronic Research Archive, 2021, 29(5): 3449-3469. doi: 10.3934/era.2021047
    [2] Messoud Efendiev, Vitali Vougalter . Linear and nonlinear non-Fredholm operators and their applications. Electronic Research Archive, 2022, 30(2): 515-534. doi: 10.3934/era.2022027
    [3] Bidi Younes, Abderrahmane Beniani, Khaled Zennir, Zayd Hajjej, Hongwei Zhang . Global solution for wave equation involving the fractional Laplacian with logarithmic nonlinearity. Electronic Research Archive, 2024, 32(9): 5268-5286. doi: 10.3934/era.2024243
    [4] Huali Wang, Ping Li . Fractional integral associated with the Schrödinger operators on variable exponent space. Electronic Research Archive, 2023, 31(11): 6833-6843. doi: 10.3934/era.2023345
    [5] Tao Zhang, Tingzhi Cheng . A priori estimates of solutions to nonlinear fractional Laplacian equation. Electronic Research Archive, 2023, 31(2): 1119-1133. doi: 10.3934/era.2023056
    [6] Hui Jian, Min Gong, Meixia Cai . Global existence, blow-up and mass concentration for the inhomogeneous nonlinear Schrödinger equation with inverse-square potential. Electronic Research Archive, 2023, 31(12): 7427-7451. doi: 10.3934/era.2023375
    [7] Xiaoming An, Shuangjie Peng . Multi-peak semiclassical bound states for Fractional Schrödinger Equations with fast decaying potentials. Electronic Research Archive, 2022, 30(2): 585-614. doi: 10.3934/era.2022031
    [8] Shasha Bian, Yitong Pei, Boling Guo . Numerical simulation of a generalized nonlinear derivative Schrödinger equation. Electronic Research Archive, 2022, 30(8): 3130-3152. doi: 10.3934/era.2022159
    [9] Ibtissam Issa, Zayd Hajjej . Stabilization for a degenerate wave equation with drift and potential term with boundary fractional derivative control. Electronic Research Archive, 2024, 32(8): 4926-4953. doi: 10.3934/era.2024227
    [10] Yang Liu, Wenke Li . A family of potential wells for a wave equation. Electronic Research Archive, 2020, 28(2): 807-820. doi: 10.3934/era.2020041
  • In this paper, we consider an HIV infection model with saturated infection rate, intracellular delay and saturated cytotoxic T lymphocyte (CTL) immune response. By calculation, we obtain immunity-inactivated reproduction number R0 and immunity-activated reproduction number R1. By analyzing the distribution of roots of the corresponding characteristic equations, we study the local stability of an infection-free equilibrium, an immunity-inactivated equilibrium and an immunity-activated equilibrium of the model. By constructing suitable Lyapunov functionals and using LaSalle's invariance principle, we show that if R0<1, the infection-free equilibrium is globally asymptotically stable; If R1<1<R0, the immunity-inactivated equilibrium is globally asymptotically stable; If R1>1, the immunity-activated equilibrium is globally asymptotically stable. Sensitivity analyses are carried out to show the effects of parameters on the immunity-activated reproduction number R1 and the viral load.


    In this paper, we examine the following Schrödinger equation:

    {iψt=(Δ)sψ+|x|2ψ|ψ|2σψinRN×[0,),ψ(0,x)=ψ0(x)Hs(RN), (1)

    where 0<s<1, σ>0, N1 and ψ:RN×[0,)C is the wave function with initial condition ψ0(x) belonging to the following Sobolev space:

    Hs(RN):={uL2(RN):RNRN|u(x)u(y)|2|xy|N+2sdydx<}.

    with

    uHs(RN)=RNRN|u(x)u(y)|2|xy|N+2sdydx+RN|u(x)|2dx.

    The fractional Laplacian (Δ)s is defined via a pseudo-differential operator

    (Δ)su(x)=F1[|ξ|2s[u]], s>0. (2)

    For the Cauchy problem, formally we have two conserved quantities [17] by multiplying the conjugate of ψ to both sides and taking the integral: The mass of the wave function:

    M(t)=||ψ(,t)||2:=RN|ψ(x,t)|2dxM(0) (3)

    and the total energy:

    E(t)=RN[Re(ψ(x,t)(Δ)sψ)+|x|2|ψ|21σ+1|ψ(x,t)|2(σ+1)]dxE(0). (4)

    In recent years, a great attention has been focused on the study of problems involving the fractional Laplacian, which naturally appears in obstacle problems, phase transition, conservation laws, financial market. Nonlinear fractional Schrödinger equations (FNLS) have been proposed by [20,19,9] in order to expand the Feynman path integral, from the Brownian like to the Lévy like quantum mechanical paths. The stationary solutions of fractional nonlinear Schrödinger equations have also been intensively studied due to their huge importance in nonlinear optics and quantum mechanics [20,19,14,12]. The most interesting solutions have the special form:

    ψ(x,t)=eiλtu(x),λR,u(x) isasmoothfunction. (5)

    They are called the standing waves. These solutions reduce (1) to a semilinear elliptic equation. In fact, after plugging (5) into (1), we need to solve the following equation

    (Δ)su(x)+|x|2u(x)|u(x)|2σu(x)=λu(x)inRN. (6)

    The case s=1 has been intensively studied by many authors (See [24]). There also exist a considerable amount of results concerning the standing waves and long time behavior of fractional nonlinear Schrödinger equations, we refer the readers to [4,7,6,10,13,16,21,23,22,1] and the references therein.

    In this paper, we mainly focus on the solutions to (6). To the best of our knowledge, our results are new and will open the way to solve other classes of fractional Schrödinger equations. This paper has two main parts: In the first part, we address the existence of standing waves through a particular variational form, whose solutions are called ground state solutions. We prove the existence of ground state solutions (Theorem 2.1), and show some qualitative properties like monotonicity and radiality (Lemma 4). We also show that the ground state solutions are orbitally stable (Def 4.2, Theorem 2.2) if we have the uniqueness of the solutions for the Cauchy problem (1) (Theorem 4.3). We have also addressed the critical case σ=2sN, which is consistent with the case s=1 in [11]. The second part of this article deals with the numerical method to solve (1) and the existence of ground state solutions as well as the optimality of our conditions. In this part, we were not only able to show the existence of ground state solutions for 0<σ<2sN but we also gave a characterizing variational problem ((65)-(66)), which was crucial to find the standing waves for the subcritical 2sNσ<2sN2s. The numerical results provided a good explanation of the effect of s on the ground state solution. To reach this goal, we showed that the ground state solution is continuous and decreasing with respect to s in L2 and L norm (Figure 3), which is a similar phenomenon to [18]. Besides, like Gross–Pitaevskii equation [11], we examined the convergence property of λc. It turned out that this convergence property also holds in our case. Second, we checked the stability of ground state solutions for different s. If we add a small perturbation to the initial condition, for different s, the absolute value of the solution will always have periodic behavior, which shows the orbital stability (Figure 7). Furthermore, surprisingly, when s becomes smaller, the stability is worse, which means the oscillation amplitude in the periodic phenomenon becomes larger (Figure 9, 10a). We then address the case where the harmonic potential is not radial, and we obtained nonradial symmetrical ground state solution (Figure 6a, 6b). Finally, we provided interesting numerical results for the time dynamics of the fractional nonlinear Schrödinger equation (FNLS).

    Figure 3.  Ground state solutions with σ=1 and different s.
    Figure 6.  Ground state solutions with non-symmetric potential.
    Figure 7.  Energy and stability check with s=0.8, σ=1.
    Figure 9.  Abosolute value of solution with different s.
    Figure 10.  Stability test with ψs0(x)=0.9us0(x).

    The main difficulty of constructing ground state solutions comes from the lack of compactness of the Sobolev embeddings for the unbounded domain RN. However, by defining an appropriate function space, in which the norm of the potential is involved, we "recuperate" the compactness (see Lemma 3.1). This fact, combined with rearrangement inequalities is the key point to prove the existence of ground state solutions. In the numerical part, the presence of the harmonic potential term is challenging. One can't take Fourier transform directly on both sides of the equation like [18] because we have the nonlinear term. Unlike in [11], we also can't use finite-difference directly since fractional Laplacian is not a local term. Consequently, we use the splitting method and adopt an idea from [8], where they numerically solve fractional Schrödinger equation without potential energy. By our splitting, we can obtain specific solutions in each small step and also preserve the mass (3). For the ground state solutions, the classical Newton's method [11] is too slow because we have to deal with fractional Laplacian. To overcome this, we borrow an idea from [2] and use normalized gradient flow (NGF) to find the ground state solutions. Moreover, for the case 2sNσ2sN2s, we have noticed that the energy in the original variational problem can not be bounded from below, therefore, we present a new characterizing variational problem ((65)-(66)) to establish the existence of ground state solutions.

    The paper is organized as follows. In section 2, we give our main results about the existence of ground state solutions and the orbital stability of standing waves. In section 3, we provide proof of existence. Then, in section 4, we discuss the orbital stability. In section 5, we use Split-Step Fourier Spectral method to solve (1) numerically. In section 6, instead of using common iterative Newton's method, we use the NGF method to find ground states when 0<σ2sN2s. Finally, in section 7, we present our numerical results for the dynamics (1) and compare them with other kinds of Schrödinger equations ([11,18]).

    We use a variational formulation to examine the solution to (6). First, note that if λ=0, we can find solutions u(x) to (6) from the critical points of the functional J:Hs(RN)R defined as:

    J(u)=12su2L2(RN)+12RN|x|2|u|2dx12σ+2RN|u|2σ+2dx, (7)

    where .L2(RN) is the L2-norm and suL2(RN) is defined by

    su2L2(RN)=CN,sRNRN|u(x)u(y)|2|xy|N+2sdxdy,

    with some normalization constant CN,s.

    We can derive (7) by multiplying smooth enough test function v(x) on both sides of (6) and taking the integral over x. However, instead of directly finding the critical points of (7), we consider a reconstructed variational problem, which can help us find solutions with different λ and any energy. Specifically, for a fixed number c>0, we need to solve the following constrained minimization problem.

    Ic=inf{J(u):uSc}, (8)

    with

    Sc={uΣs(RN):RN|u|2dx=c2}, (9)

    where

    Σs(RN)={uHs(RN):uΣs(RN):=u2L2(RN)+su2L2(RN)+xu2L2(RN)<}. (10)

    is a Hilbert space, with corresponding natural inner product.

    We claim that for each minimizer u(x) of the constrained minimization problem (8), there exists some λ such that (u(x),λ) is a solution to (6). To prove the claim, we consider λ as a Lagrange multiplier and define

    J(u)=J(u)+λ(u2L2(RN)c2). (11)

    The minimizer to problem (8) must be the critical point of (11), satisfying:

    J(u)u=0 (12)

    and

    J(u)λ=0, (13)

    where (12) implies (6) and (13) implies (9). In this paper, we will mainly focus on the minimizers of problem (8). The following theorem discusses the existence of such minimizers.

    Theorem 2.1. If 0<σ<2sN, then (8) admits a nonnegative, radial and radially decreasing minimizer.

    Remark 1. There were some works considering standing waves for fractional nonlinear Schrödinger equation with different potential enery and nonlinearity terms. In [6], the authors focus on the case where the nonlinearity term is asymptotically linear with ψ. In [7], the author considers strictly positive potential energy. In our case, the nonlinearity term increases polynomially with ψ and potential energy is not strictly positive, therefore previous technics can not be applied here directly.

    Remark 2. The condition 0<σ<2sN is important in our proof for the existence of minimizers. For the critical case σ=2sN, we were able to obtain interesting results (section 7).

    After we construct the ground state solutions, we further investigate their stability. By the definition of (5), the ground state solution moves around a circle when the time changes. Therefore we consider and prove the orbital stability of ground state solution (Def 4.2).

    Theorem 2.2. Suppose that 0<σ<2sN and (1) has a unique solution with conserved mass (3) and energy (4), then the ground state solutions constructed in Theorem 2.1 are orbitally stable.

    In this section, we will establish the existence of ground state solutions of (6), the main difficulty comes from the lack of compactness of the Sobolev embeddings. Usually, at least when the potential in (1) is radially symmetric and radially increasing, such difficulty is overcome by considering the appropriate function space. More precisely, we have

    Lemma 3.1. Let 2p<2NN2s, then the embedding Σs(RN)Lp(RN) is compact.

    Proof. First, when p=2, choose R>0, then for any u(x)0Σs(RN), we have

    |x|R|u|2dx|R|2|x|R|x|2|u|2dx<|R|2u2Σs(RN) (14)

    By the classical Sobolev embedding theorem, for any fixed R, Hs(|x|<R) is compactly embedded in L2(|x|<R). Therefore, for any bounded sequence in Σs, we choose Rn>0. And for each n, we pick out the subsequence that converges in L2(|x|<Rn) from former convergence sequence in L2(|x|<Rn1). Finally, using the diagonal method combined with (14), we find the convergence sequence in L2(RN).

    Second, for p>2, using the fractional Gagliardo-Nirenberg inequality [14], we have

    upLp(RN)Ku1θL2(RN)suθL2(RN)Ku1θL2(RN)uθΣs(RN), (15)

    for some positive constant K, where θ=N(p2)2sp<1. Since Σs(RN) is compactly embedded into L2(RN), (15) directly implies Σs(RN) is compactly embedded into Lp(RN) for 2p<2NN2s.

    Then we have a lemma showing the existence of Ic and the boundedness of the minimizing sequences.

    Lemma 3.2. If 0<σ<2sN, then Ic> and all the minimizing sequences of (8) are bounded in Σs(RN).

    Proof. First, we prove that J(u) is bounded from below. Using the fractional Gagliardo-Nirenberg inequality [14], we certainly have

    uL2σ+2(RN)Ku1θL2(RN)suθL2(RN), (16)

    for some positive constant K, where θ=Nσ2s(σ+1).

    On the other hand, let ϵ>0, and p,q>1 such that 1p+1q=1, then, using Young's inequality, one gets that

    u(2σ+2)(1θ)L2(RN)suθ(2σ+2)L2(RN)1pϵpsupθ(2σ+2)L2(RN)+1qϵquq(1θ)(2σ+2)L2(RN). (17)

    Combining (16) and (17), we obtain that for any uSc,

    RN|u(x)|2σ+2dxϵpK2σ+2psu2L2(RN)+K2σ+2qϵqc2q(1θ)(1+σ), (18)

    where p=1θ(1+σ)=2sNσ, q=11θ(1+σ).

    Hence, from (18) we get:

    J(u) (19)
    =12su2L2(RN)+12RN|x|2|u|2dx12σ+2RN|u|2σ+2dx (20)
    12su2L2(RN)+12RN|x|2|u|2dx (21)
    12σ+2(ϵpK2σ+2psu2L2(RN)+K2σ+2qϵqc2q(1θ)(1+σ)) (22)
    (12ϵpK2σ+22p(σ+1))su2L2(RN)+12RN|x|2|u|2dxK2σ+2c2q(1θ)(1+σ)2q(1+σ)ϵq. (23)

    Then we choose ϵ small enough in (23) to make (12ϵpK2σ+22p(σ+1))>0, which implies that Ic>. Finally, for all minimizing sequences {un}, since J(un) is bounded from below, {un} is bounded in Σs(RN) by (23).

    Now, we can use compactness (Lemma 3.1) and boundedness (Lemma 3.2) to prove our existence Theorem 2.1.

    Proof. Let {un} be a minimizing sequence of (8). By Lemma 3.2, {un} is bounded in Σs(RN). Up to a subsequence, there exists u such that un converges weakly to u in Σs(RN).

    Since 2σ+2<2NN2s and Σs(RN) is compactly embedded in Lp(RN) for any p such that 2p<2NN2s, we can further prove that un will converge strongly to u in L2(RN) and L2σ+2(RN) (Lemma 3.1). In particular, unu in L2(RN), which implies that uSc.

    On the other hand, thanks to the lower semi-continuity, we have

    xuL2(RN)+suL2(RN)lim infnsunL2(RN)+xunL2(RN).

    Therefore

    IcJ(u)lim infnJ(un)=Ic, (24)

    which yields that u is a minimizer.

    The second step consists of constructing a nonnegative, radial, and radially decreasing minimizer. First, we note that:

    s|u|L2(RN)suL2(RN), (25)

    which implies that J(|u|)J(u). Using the Schwarz symmetrization [15], we construct a symmetrization function u, which is a radially-decreasing function from RN into R with the property that

    meas{xRN:|u(x)|>μ}=meas{xRN:u(x)>μ}foranyμ>0.

    It is well-known [15] that

    {RN|u|2σ+2dx=RN|u|2σ+2dx,RN|u|2dx=RN|u|2dx,RN|x|2|u|2dxRN|x|2|u|2dx. (26)

    Besides, from [3] Lemma 6.3, we also have

    suL2(RN)suL2(RN). (27)

    Combining (26) and (27), we obtain that

    J(|u|)J(|u|)J(u),forany uΣs(RN).

    Remark 3. By (24) and the weakly convergence, we can also see that:

    xuL2(RN)+uL2(RN)+suL2(RN)=limnsunL2(RN)+unL2(RN)+xunL2(RN),

    which implies that there is a minimizing subseqence unk converging to u in Σs(RN).

    Remark 4. If uΣs(RN) is a minimizer to (8), we must have

    Ic=J(u)=J(|u|)=J(|u|). (28)

    By (27) and (28), we obtain that

    suL2(RN)=s|u|L2(RN), (29)

    which implies that u doesn't change sign almost everywhere. Furthermore, we also have

    RN|x|2(|u|)2dx=RN|x|2|u|2dx. (30)

    which further implies |u|=|u| a.e. by [11]. In conlusion, we obtain a minimizer u that doesn't change sign and is radially decreasing almost everywhere.

    In this section, we will deal with the orbital stability of the ground state solutions. Let us introduce the appropriate Hilbert space:

    ˜Σs(RN):={ω=u+iv:(u,v)Σs(RN)×Σs(RN)},

    equipped with the norm ω2˜Σs(RN)=u2Σs(RN)+v2Σs(RN), which is a Hilbert space.

    In term of the new coordinates, the energy functional reads

    ˜J(ω)=12sω2L2(RN)+12RN|x|2|ω(x)|2dx12σ+2RN|ω|2σ+2dx,

    where sω2L2(RN)=su2L2(RN)+sv2L2(RN), we can also prove that ˜J(ω) remains as a constant with time t if ω(t,x) is a solution to (1).

    Then, for all c>0, we set a similar constrained minimization problem

    ˜Ic=inf{˜J(ω),ω˜Sc},

    where ˜Sc is defined by:

    ˜Sc={ω˜Σs(RN),RN|ω(x)|2dx=c2}.

    We also introduce the following sets

    Oc={uSc:Ic=J(u)},˜Oc={ω˜Sc:˜Ic=˜J(ω)}.

    Proceeding as in [4,16], we have the following lemma:

    Lemma 4.1. If 0<σ<2sN, then the following properties hold true:

    (i) The energy functional J and ˜J are of class C1 on Σs(RN) and ˜Σs(RN) respectively.

    (ii) There exists a constant C>0 such that

    J(u)Σ1s(RN)C(uΣs(RN)+u2σ+1Σs(RN)),˜J(ω)Σ1s(RN)C(ωΣs(RN)+ω2σ+1Σs(RN)).

    (iii) All minimizing sequences for ˜Ic are bounded in ˜Σs(RN) and all minimizing sequences for Ic are bounded in Σs(RN).

    (iv) The mappings cIc,˜Ic are continuous.

    (v) Any minimizing sequence of Ic, ˜Ic are relatively compact in Σs(RN), ~Σs(RN).

    (vi) For any c>0,

    Ic=˜Ic.

    Proof. (i) We follow the steps of Proposition 2.3 [16] by choosing g(x,t)=tσ. For any u,vΣs(RN), we can see that the last term of functional

    RN|u(x)|2σu(x)v(x)dx

    is of class C1 on Σs(RN). Then by the definition of Σs(RN) (see (10)), the first two terms of the functional are of class C1 on Σs(RN).

    (ii) From (i), J is of class C1 on Σs(RN). Moreover, for all u,vΣs(RN), we have

    J(u),v=CN,SRNRN|u(x)u(y)||v(x)v(y)||xy|N+2sdxdy+RN|x|2u(x)v(x)dxRN|u(x)|2σu(x)v(x)dx.

    For the last term, by Hölder's inequality, we have

    RN|u(x)|2σu(x)v(x)dxu2σ+1L2σ+2(RN)vL2σ+2(RN).

    Therefore, there exists C>0 such that

    J(u)Σ1s(RN)C(uΣs(RN)+u2σ+1Σs(RN)).

    (iii) This is a direct result of Lemma 3.2.

    (iv) Let c>0 and {cn}(0,) such that cnc. It suffices to prove that IcnIc. By the definition of Icn, for any n, there exists unScn such that

    IcnJ(un)<Icn+1n. (31)

    From (iii), there exists a constant C1>0 such that for all n, we have

    unΣs(RN)C1,nN.

    Setting vn=ccnun, then we have that for all nN

    vnScandunvnΣs(RN)=|1ccn|unΣs(RN)C1|1ccn|,

    which implies that

    unvnΣs(RN)C1+1fornlargeenough. (32)

    We deduce by part (ii) that there exists a positive constant K:=K(C1) such that

    J(u)Σ1s(RN)K,foralluΣs(RN)withuΣs(RN)2C1+1. (33)

    By (32) and (33), we get that

    |J(vn)J(un)|=|10ddtJ(tvn+(1t)un)dt|supuΣs(RN)2C1+1J(u)Σ1s(RN)vnunΣs(RN)KC1|1ccn|. (34)

    Then, from (31) and (34), we obtain that

    IcnJ(un)1nJ(vn)KC1|1ccn|1nIcKC1|1ccn|1n

    Combining this with the fact that limncn=c, it yields

    lim infnIcnIc. (35)

    Now, from Lemma 3.2 and by the definition of Ic, there exists a positive constant C2 and a sequence {un}Sc such that

    unΣs(RN)C2andlimnJ(un)=Ic.

    Set vn=cncun, then vnScn, there exists a constant L=L(C2) such that

    vnunΣs(RN)C2|1cnc|and|J(vn)J(un)|LC2|1cnc|.

    Combining this with (31), we obtain that

    IcnJ(vn)J(un)+LC2|1cnc|.

    Since limncn=c, we have

    lim supnIcnIc. (36)

    It follows from (35) and (36) that

    limnIcn=Ic.

    (v) This is a direct result of Remark 3.

    (vi) First, we can see Σs(RN)˜Σs(RN). And for any ωΣs(RN), we have

    ˜J(ω)=J(ω),

    which implies that

    Ic˜Ic. (37)

    Second, for any ω˜Σs(RN), we have

    sω2L2(RN)s|ω|2L2(RN),

    which implies that

    ˜J(ω)˜J(|ω|)=J(ω)Ic,ω˜Σs(RN),

    and

    ˜IcIc. (38)

    Combining (37) and (38), we finally have ˜Ic=Ic.

    Now, for a fixed c>0, we use the following definition of stability (see [5]):

    Definition 4.2 We say that ˜Oc is stable if

    ˜Oc is not empty.

    ● For all ω0˜Oc and ε>0, there exists δ>0 such that for all ψ0˜Σs(RN), we have

    ω0ψ0˜Σs(RN)<δinfω˜Ocωψ˜Σs(RN)<ε,

    where ψ denotes the solution of (1) corresponding to the initial data ψ0.

    If ˜Oc is stable, we say that the ground state solutions in ˜Oc are orbitally stable. The following theorem states the orbital stability of ˜Oc.

    Theorem 4.3. Suppose that 0<σ<2sN, and (1) with initial data ψ0˜Σs(RN) has a unique solution ψ(t,x)˜Σs(RN) with

    ψ(t,.)L2(RN)=ψ0(t,.)L2(RN)and˜J(ψ(t,.))=˜J(ψ0(t,.)), (39)

    then ˜Oc is stable.

    Proof. The proof is by contradiction: Suppose that ˜Oc is not stable, then there exists ϵ0>0,ω0˜Oc and a sequence Φn0˜Σs(RN) such that ω0Φn0˜Σs(RN)0asn, but

    infZ˜OcΦn(tn,.)Z˜Σs(RN)ε (40)

    for some sequence {tn}R, where Φn(t,.) is the unique solution of problem (1) corresponding to the initial condition Φn0.

    Let ωn=Φn(tn,.)=(un,vn)˜Σs(RN). Since ω˜Sc and ˜J(ω)=˜Ic, it follows from the continuity of .L2(RN) and ˜J in ˜Σs(RN) that

    Φn0L2(RN)cand˜J(Φn0)˜Ic, n.

    Thus, we deduce from (39) that

    ωnL2(RN)=Φn0L2(RN)cand˜J(ωn)=˜J(Φn0)˜Ic, n. (41)

    Since {ωn}˜Σs(RN), it is easy to see that {|ωn|}Σs(RN). On the other hand, Lemma 4.1 (iii) and the proof of Lemma 3.2 imply that {ωn} is bounded in ˜Σs(RN). Hence, by passing to a subsequence, there exists ω=(u,v)˜Σs(RN) such that

    unu,vnvandlim infnsunL2(RN)+svnL2(RN)exists. (42)

    Now, by a straightforward computation, we obtain that

    ˜J(ωn)˜J(|ωn|)=12sωn2L2(RN)12s|ωn|2L2(RN)0, (43)

    which implies that

    ˜Ic=limn˜J(ωn)lim supnJ(|ωn|).

    Besides, by (41),

    ωn2L2(RN)=|ωn|2L2(RN)=c2nc2.

    It follows from Lemma 4.1 that we have

    lim infnJ(|ωn|)lim infnIcn=Ic.

    Hence

    ˜Ic=limn˜J(ωn)=limnJ(|ωn|)=Ic. (44)

    It follows from (42), (43) and (44) that

    limnsun2L2(RN)+svn2L2(RN)s(u2n+v2n)122L2(RN)=0,

    which is equivalent to say that

    limnswn2L2(RN)=limns|wn|2L2(RN). (45)

    The boundedness of ωn in ~Σs(RN) and (45) imply that |wn| is bounded in Σs(RN). Using a similar argument to Lemma 3.1, there exists φΣs(RN) such that

    |ωn|φinΣs(RN)andφL2(RN)=cwithJ(φ)=Ic. (46)

    Next, let us prove φ=|ω|=(|u|2+|v|2)1/2. Using (42), it follows that

    unuandvnvinL2(B(0,R))forallR>0.

    Since |(u2n+v2n)12(u2+v2)12||unu|2+|vnv|2, one has

    (u2n+v2n)12(u2+v2)12inL2(B(0,R)).

    But |ωn|=(u2n+v2n)12φinΣs(RN)L2(RN). Thus, we certainly have

    (u2+v2)12=|ω|=φ.

    This further implies that

    ωL2(RN)=φL2(RN)=c,ωL2σ+2(RN)=φL2σ+2(RN). (47)

    and

    12RN|x|2|ω|2dx12σ+2RN|ω|2σ+2dx=12RN|x|2|φ|2dx12σ+2RN|φ|2σ+2dx=limn12RN|x|2|ωn|2dx12σ+2RN|ωn|2σ+2dx. (48)

    Additionally, by the lower semi-continuity, we further have

    ω2L2(RN)lim infnωn2L2(RN). (49)

    Combining (48) with (49) and ω˜Sc, we finally obtain that

    ˜Ic˜J(ω)limn˜J(ωn)=˜Ic,

    which implies that

    ω˜Ocandsω2L2(RN)=limnsωn2L2(RN). (50)

    Finally, by (46), (47), and (50), we obtain that

    ωnωin˜Σs(RN),

    which contradicts to (40).

    In this section, we consider numerical methods to solve (1) and introduce the Split-Step Fourier Spectral method.

    First, we truncate (1) into a finite computational domain [L,L]N with periodic boundary conditions:

    {iψt=(Δ)sψ+|x|2ψ|ψ|2σψ, t>0,ψ(0,x)=ψ0(x), (51)

    for x[L,L]N.

    Let τ>0 be the time step, and define the time sequence tn=nτ for n0 and the mesh size h=2L/J, where J is a positive even integer. The spatial grid points are

    (xj)n=L+(j)nh, 1nN, (52)

    where j is a N-dimension integer vector with each component between 0 and J.

    Denote ψnj the numerical approximation of the solution ψ(xj,tn). By the definition of fractional Laplacian in (2), we use the Fourier spectral method for spatial discretization. Hence, we assume the ansatz:

    ψ(xj,t)=kK^ψk(t)exp(iμkxj), (53)

    where K={kRN|J/2klJ/21,1lN}, (μk)k=kkπ/L, 1kN.

    Now, we introduce the Split-step Fourier Spectral method. The main idea of this method is to solve (51) in two splitting steps from t=tn to t=tn+1 :

    iψt=|x|2ψ|ψ|2σψ, (54)
    iψt=(Δ)sψ. (55)

    First, by multiplying ψ on both sides of (54) and subtracting it from its conjugate, we obtain that |ψ(x,t)|=|ψ(x,tn)| for any t[tn,tn+1). Then (54) can be simplified to

    iψt=|x|2ψ|ψ(x,tn)|2σψ. (56)

    Second, taking Fourier transform on both sides of (55), we obtain that

    id^ψk(t)dt=|μk|2s^ψk(t). (57)

    We use the second order Strang splitting method with (56) and (57) as follows:

    ψn,1j=ψnjexp(i(|xj|2|ψnj|2σ)τ/2) (58)
    ψn,2j=kK^ψn,1kexp(i|μk|2sτ)exp(iμkxj) (59)
    ψn+1j=ψn,2jexp(i(|xj|2|ψn,2j|2σ)τ/2) (60)

    where j comes from (52) and n0. For n=0, initial condition (51) is discretized as:

    ψ0j=ψ0(xj) (61)

    This method has spectral order accuracy in space and second order in time. Similar to [8], this method preserves discrete mass corresponding to (3) defined as

    Mn=(hNj|ψnj|2)1/2. (62)

    To find ground state solutions, we have to solve the following equation corresponding to u(x):

    (Δ)su+|x|2u|u|2σu=λu xRN. (63)

    As discussed previously, for σ<2sN, we can solve (8)-(10) to find a solution to (63). In order to calculate the minimizer of J(u) in Sc, we use the normalized gradient flow method (NGF) [2]. We first apply the steepest gradient decent method to the energy functional J(u) without constraint. Then we project the solution back onto the sphere Sc to make sure that the constraint ||u||L2(RN)=c is satisfied.

    Thus, for a given sequence of time 0=t0<t1<...<tn with fixed time step τ, we compute the approximated solution u(n) of the partial differential equation

    ut=E(u)u

    and project it onto Sc at each step. Specifically,

    {˜ut=(Δ)s˜u|x|2˜u+|˜u|2σ˜u,tn<t<tn+1,˜u(x,tn)=u(n)(x),u(n+1)(x)=c˜u(x,tn+1)||˜u(,tn+1)||L2(RN),.

    Here, we use semi-implicity time discretization scheme:

    {˜u(n+1)˜u(n)Δt=(Δ)s˜u(n+1)|x|2˜u(n+1)+|˜u(n)|2σ˜u(n+1),tn<t<tn+1,˜u(x,tn)=u(n)(x),u(n+1)(x)=c˜u(n+1)(x)||˜u(n+1)()||L2(RN),

    with

    δsxu|j=kK|μk|2sexp(iμkxj)

    to discretize fractional Laplacian, where K, μk are defined in (53).

    In conlusion, in each step, we solve:

    {˜u(n+1)ju(n)jτ=δsx˜u(n+1)|j(|xj|2|u(n)j|2σ)˜u(n+1)ju(n+1)j=c˜u(n+1)jMn, (64)

    where Mn is defined in (62), j comes from (52) and n0. For n=0, we guess a starting function and discretize it as (61).

    We need to notice that we can only solve (8) for σ<2sN. If σ2sN, uL2σ+2(RN) can not be bounded by Σs(RN), which causes Ic= in Sc. However, we can use another characterizing variational form and find standing waves to (1) for 2sNσ<2sN2s.

    For σ<2sN2s, we define the following constrained minimization problem:

    Lc=inf{K(u):uTc} (65)

    with

    Tc={uΣs(RN):uL2σ+2(RN)=c},
    K(u)=12su2L2(RN)+12RN|x|2|u(x)|2dx+ω2RN|u(x)|2dx, (66)

    where Σs(RN) is defined as (10).

    For uTc, we have the estimate

    RN|u(x)|2dx|R|2|x|>R|x|2|u|2dx+|x|<R|u|2dx|R|2|x|>R|x|2|u|2dx+CR,σu2L2σ+2(RN),

    where CR,δ depends on R,δ. Hence, for ω<0, if we choose R large enough to make (12+ωR2)>0, we obtain that

    K(u)12su2L2+(12+ωR2)RN|x|2|u(x)|2dx+ω2CR,σc2> (67)

    for uTc. Besides, if ω>0, (67) is greater than 0. Similar to Lemma 3.2 and Theorem 2.1, there exists a local minimizer for (65) with any ω and c.

    Now we see λ as the Lagrange multiplier like (11) but with a different functional

    K(u)=K(u)λ2σ+2(u2σ+2L2σ+2(RN)c2σ+2),

    then the critical points u and λ satisfy

    (Δ)su+|x|2uλ|u|2σu+ωu=0,xRN (68)

    by K(u)u=0. If ω>0, by multiplying ¯u on both sides of (68) and taking the integral, we can see that:

    K(u)=K(u)=λ2c2σ+2>0, (69)

    which implies that λ>0. Defining uω,c=λ1/2σu, when ω>0, we obtain that

    (Δ)suω,c+|x|2uω,c|uω,c|2σuω,c+ωuω,c=0,xRN, (70)

    which means uω,ceiωt is one standing wave solution to (1). We need to mention that (69), (70) actually showed that we can find a ground state solution by solving (66) if Lc=K(u)>0. In fact, we have K(u)>0 with ω<0 but not very small. This is related with the smallest eigenvalue of (Δ)s+|x|2 (see detail [11]).

    Now, for 2sNσ<2sN2s, we use the NGF method and the semi-implicity time discretization scheme to solve the constrained problem (65). Similar to (64), the scheme is

    {˜u(n+1)ju(n)jτ=δsx˜u(n+1)|j(|xj|2+u(n)j)˜u(n+1)ju(n+1)j=c˜u(n+1)jMn2σ+2, (71)

    where Mn2σ+2 is discrete L2σ+2(RN) norm

    Mn2σ+2=(hNj|unj|2σ+2)1/(2σ+2)

    and j comes from (52) for n0. For n=0, we guess a starting function and discretize it as (61).

    In this section, we show some numerical results, which can help us understand the ground state solution and also illustrate theoretical results. We have mainly investigated: 1. Ground state solutions with different s. 2. Ground state solutions with non-symmetric potentials. 3. Stability and dynamics.

    First, we solve (8) numerically by (64) in one dimension N=1 for the case s=0.8 and σ=1 to obtain a ground state solution u0(x). From figure 1a we see that u0(x) is radially decreasing as Remark 4.

    Figure 1.  Ground state solution with s=0.8, σ=1, L=10 and J=5000.

    Second, we use u0(x) as initial condition in (1) and investigate time evolution of the standing wave (figure 1b, 2a, 2b). As expected, we see that |ψ(x,t)| is conserved. And the real and imaginary part of the solution change periodicly with time t.

    Figure 2.  Time dynamics of standing waves with s=0.8, σ=1, L=10 and J=5000.

    By Theorem 2.1, we can obtain the existence of ground state solutions with σ<2sN. We change s but keep σ=1 to obtain ground state solutions with different s. From figure 3a, we can see when s approaches to 0.5, ground state solutions become peaked with faster spatial decay. This is a similar result to the case without potential [18]. We also check Σ1(R) of ground state solutions when s0.5 in (figure 3b), whose growth shows regularity of the solution becomes worse.

    From figure 3a and figure 3b, it seems ground state solutions change continuously with s. We use L2 distance between us0(x)u1(x) to check and see the convergence of ground state solutions in L^2 space with s\rightarrow 1 (figure 4).

    Figure 4.  L^2 distance between ground state solutions of s<1 and s = 1 with \sigma = 1 .

    Then, we test another two things. The first is the relationship between the constrained minimal energy in (8) and s . We calculate the discrete energy by

    \begin{equation} E(s) = h{\sum^{J-1}\limits_{j = 0}}\left[{\sum^{J/2-1}\limits_{l = -J/2}}|\mu_l|^{2s}|\widehat{u}_l|^2+|x_j|^2|u_j|^2-\frac{1}{\sigma+1}|u_j|^{2(\sigma+1)}\right]. \end{equation} (72)

    From figure 5a, we find the energy's dependence (E(s)) on s is monotonic. When s approaches to 0.5 , the energy will approach to -\infty because of the nonlinear term. There are two reasons. First, we keep the L^2 norm of u (we test with the same mass c ), but the potential term becomes small since u gathers around 0 . Second, in Lemma 3.2, we need \sigma<\frac{2s}{N} to bound \|u\|_{L^{2\sigma+2}(\mathbb{R}^N)} , whose boundedness becomes worse when s\rightarrow0.5 . This is different from [18], where they didn't use the variation form and keep the L^2 norm.

    Figure 5.  Energy and \lambda_c .

    Second, we test the relationship between c and \lambda_c , where c , \lambda_c are mass and Lagrange multiplier corresponding to the minimizer (8)-(10). By [11], in the case s = 1 , there exists \lambda_0 such that for any \sigma ,

    \lim\limits_{c\rightarrow 0}\lambda_c = \lambda_0\, .

    We also test this with s = 0.8 . From figure 5b, for different \sigma , when c\rightarrow 0 , \lambda_c will also converge to a same value \lambda_0 .

    Up to now, we only consider the case with radially symmetrical potential. However, when the potential is not radially symmetric, we can still find a standing wave to (1) using (8). We tried cases where the potential has the form |x|^2+a\sin(2\pi x) with a = 1 and a = 5 . From figure 6a and 6b, if we add a nonsymmetrical perturbation to potential, we won't get radially symmetrical ground state solutions.

    First, we consider the case s>\frac{N\sigma}{2} , which is covered by our Theorem 2.2. From the theoretical results and figure 1b, we see that standing waves preserve |\psi(x, t)| with time t . Therefore, we use |\psi(x, t)| to draw graphs and test its stability. We consider the case s = 0.8 and \sigma = 1 . We first test condition (39) in Theorem 2.2 with the initial condition \psi_0(x) = (1+e)*u_0(x) . Because the scheme is mass preserving, it suffices to test energy preservation, which is showed in figure 7a, 7b. Then, we test the stability of the solution, where e is a constant number. From figure 7c, 7d, when e = 0.05 , the solution almost preserves |\psi(x, t)| as we desired. When e = 0.2 , the solution shows large perturbation but still has the periodic behavior, similar to [11] and [18]. In figure 8, we compare the \|\cdot\|_{\widetilde{\Sigma}_s(\mathbb{R})} distance between ground state solutions and perturbed solutions. We can see from figure 3b and theoretical results that when s\rightarrow \frac{\sigma N}{2} , the regularity of ground state solutions becomes worse. This inspires us to investigate its stability relationship with s . By Theorem 2.2, Def 4.2, the orbital stability means we can find \omega\in\widetilde{\mathcal{O}}_{c} such that \|\omega-\psi\|_{\widetilde{\Sigma}_{s}(\mathbb{R})} is small when we only have a small perturbation in initial condition. This definition is hard to measure. Instead of checking the exact definition of the orbital stability, we test classical stability by comparing the distance between perturbed solutions and ground state solutions using the normalized \|\cdot\|_{\widetilde{\Sigma}_s(\mathbb{R})} distance:

    D(s, t) = \frac{\|u^s(x, t)-\psi^s(x, t)\|_{\widetilde{\Sigma}_{s}(\mathbb{R})}}{\|u^s(x, t)\|_{\widetilde{\Sigma}_{s}(\mathbb{R})}}.
    Figure 8.  \|\psi^s(x, t)-u^s(x, t)\|_{\widetilde{\Sigma}_s(\mathbb{R})} when s = 0.8 and \sigma = 1 .

    We test initial conditions \psi_0(x) = 0.9*u_0(x) with s = 0.8 and s = 0.6 . From figure 9, 10a, as expected, when s is small, its stability seems worse. To be complete, we test D(s, 1) with s between 0.51 and 1 in figure 10b. We can see D(1, s) increases when s approaches to 0.5 , which implies worse stability.

    Second, we try to obtain some numerical results when we touch the critical point s = \frac{\sigma N}{2} . In this case, we can't find the ground state solution through (8) because I_c = -\infty . However, as we discussed before, we can find a ground state solution related to another constrained minimization problem (65). Here, we try to use the NGF method to find a ground state solution with s = 0.5 , \sigma = 1 . We first tried positive \omega , but the projection step dominated the process (71). Therefore, we tried \omega = -0.5 and find the method does converge to a solution. From figure 11b-11d, we can see |\psi(x, t)| is almost preserved with time t and it has periodical real and imaginary part. We use it to test the finite blowup phenomenon ( \psi_0(x) = 2u_0(x) ) appeared in the case without potential [18], and this also happens with potential (figure 11e). We note here that we still can't find a perfect ground state solution, which is caused by the stability of (1) is very bad when s = 0.5 .

    Figure 11.  Standing wave, ground state solution and blow up solution with s = 0.5 , \delta = 1 .

    Finally, we test some simple time dynamics of FNLS, we let \psi_0(x) = u_0(x)e^{ikx} , which changes its phase but not absolute value. If s = 0.8 and k = 1, 20 (Figure 12a, 12b), the maximum point of |\psi(x, t)| will move along x periodically. We also test the L^\infty norm of \psi(x, t) (Figure 13a). We find that it decreases first and then approaches to some limits, which is similar to the case without potential [18].

    Figure 12.  Dynamics of FNLS.
    Figure 13.  Evolution of L^\infty .


    [1] A. M. Elaiw, Global properties of a class of HIV models, Nonlinear Anal. RWA, 11 (2010), 2253-2263. doi: 10.1016/j.nonrwa.2009.07.001
    [2] L. Wang, M. Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T cells, Math. Biosci., 200 (2006), 44-57. doi: 10.1016/j.mbs.2005.12.026
    [3] L. Cai, X. Li, M. Ghosh, B. Guo, Stability analysis of an HIV/AIDS epidemic model with treatment, J. Comput. Appl. Math., 229 (2009), 313-323. doi: 10.1016/j.cam.2008.10.067
    [4] J. Xu, Y. Geng, Y. Zhou, Global dynamics for an age-structured HIV virus infection model with cellular infection and antiretroviral therapy, Appl. Math. Comput., 305 (2017), 62-83.
    [5] M. A. Nowak, S. Bonhoeffer, G. M. Shaw, R. M. May, Anti-viral drug treatment: dynamics of resistance in free virus and infected cell populations, J. Theor. Biol., 184 (1997), 203-217. doi: 10.1006/jtbi.1996.0307
    [6] M. Y. Li, L. Wang, Backward bifurcation in a mathematical model for HIV infection in vivo with anti-retroviral treatment, Nonlinear Anal. RWA, 17 (2014), 147-160. doi: 10.1016/j.nonrwa.2013.11.002
    [7] G. P. Samanta, Permanence and extinction of a nonautonomous HIV/AIDS epidemic model with distributed time delay, Nonlinear Anal. RWA, 12 (2011), 1163-1177. doi: 10.1016/j.nonrwa.2010.09.010
    [8] M. A. Nowak, C. R. M. Bangham, Population dynamics of immune responses to persistent viruses, Science, 272 (1996), 74-79. doi: 10.1126/science.272.5258.74
    [9] R. De Boer, Which of our modeling predictions are robust?, PLoS Comput. Biol., 8 (2012), e1002593.
    [10] C. Jiang, W. Wang, Complete classification of global dynamics of a virus model with immune responses, Discrete Contin. Dyn. Syst. Ser. B, 19 (2014), 1087-1103.
    [11] A. V. M. Herz, S. Bonhoeffer, R. M. Anderson, R. M. May, M. A. Nowak, Viral dynamics in vivo: Limitations on estimates of intracellular delay and virus decay, Proc. Natl. Acad. Sci. USA, 93 (1996), 7247-7251. doi: 10.1073/pnas.93.14.7247
    [12] A. M. Elaiw, I. A. Hassanien, S. A. Azoz, Global stability of HIV infection models with intracellular delays, J. Korean Math. Soc., 49 (2012), 779-794. doi: 10.4134/JKMS.2012.49.4.779
    [13] B. Reddy, J. Yin, Quantitative intracellular kinetics of HIV type 1, AIDS Res. Hum. Retrovir., 15 (1999), 273-283. doi: 10.1089/088922299311457
    [14] A. S. Perelson, A. U. Neumann, M. Markowitz, J. M. Leonard, D. D. Ho, HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time, Science, 271 (1996), 1582-1586.
    [15] J. E. Mittler, B. Sulzer, A. U. Neumann, A. S. Perelson, Influence of delayed viral production on viral dynamics in HIV-1 infected patients, Math. Biosci., 152 (1998), 143-163. doi: 10.1016/S0025-5564(98)10027-5
    [16] D. Ebert, C. D. Zschokke-Rohringer, H. J. Carius, Dose effects and density-dependent regulation of two microparasites of Daphnia magna, Oecologia, 122 (2000), 200-209. doi: 10.1007/PL00008847
    [17] R. R. Regoes, D. Ebert, S. Bonhoeffer, Dose-dependent infection rates of parasites produce the Allee effect in epidemiology, Proc. R. Soc. Lond. Ser. B, 269 (2002), 271-279. doi: 10.1098/rspb.2001.1816
    [18] X. Song, A. U. Neumann, Global stability and periodic solution of the viral dynamics, J. Math. Anal. Appl., 329 (2007), 281-297. doi: 10.1016/j.jmaa.2006.06.064
    [19] J. K. Hale, S. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993.
    [20] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29-48. doi: 10.1016/S0025-5564(02)00108-6
    [21] M. Stafford, L. Corey, Y. Cao, E. Daar, D. Ho, A. Perelson, Modeling plasma virus concentration during primary HIV infection, J. Theor. Biol., 203 (2000), 285-301. doi: 10.1006/jtbi.2000.1076
    [22] K. Pawelek, S. Liu, F. Pahlevani, L. Rong, A model of HIV-1 infection with two time delays: Mathematical analysis and comparison with patient data, Math. Biosci., 235 (2012), 98-109. doi: 10.1016/j.mbs.2011.11.002
    [23] J. Wang, M. Guo, X. Liu, Z. Zhao, Threshold dynamics of HIV-1 virus model with cell-to-cell transmission, cell-mediated immune responses and distributed delay, Appl. Math. Comput., 291 (2016), 149-161.
    [24] S. Marino, I. B. Hogue, C. J. Ray, A methodology for performing global uncertainty and sensitivity analysis in systems biology, J. Theor. Biol., 254 (2008), 178-196. doi: 10.1016/j.jtbi.2008.04.011
  • This article has been cited by:

    1. Mona M. Saber, Nada Monir, Azza S. Awad, Marwa E. Elsherbiny, Hala F. Zaki, TLR9: A friend or a foe, 2022, 307, 00243205, 120874, 10.1016/j.lfs.2022.120874
  • Reader Comments
  • © 2021 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(5162) PDF downloads(223) Cited by(9)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog