Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article Special Issues

Proteomic and transcriptomic analyses of Cutibacterium acnes biofilms and planktonic cultures in presence of epinephrine

  • Received: 04 March 2024 Revised: 15 May 2024 Accepted: 21 May 2024 Published: 27 May 2024
  • Transcriptomic and proteomic analysis were performed on 72 h biofilms of the acneic strain Cutibacterium acnes and planktonic cultures in the presence of epinephrine. Epinephrine predominantly downregulated genes associated with various transporter proteins. No correlation was found between proteomic and transcriptomic profiles. In control samples, the expression of 51 proteins differed between planktonic cultures and biofilms. Addition of 5 nM epinephrine reduced this number, and in the presence of 5 µM epinephrine, the difference in proteomic profiles between planktonic cultures and biofilms disappeared. According to the proteomic profiling, epinephrine itself was more effective in the case of C. acnes biofilms and potentially affected the tricarboxylic acid cycle (as well as alpha-ketoglutarate decarboxylase Kgd), biotin synthesis, cell division, and transport of different compounds in C. acnes cells. These findings are consistent with recent research on Micrococcus luteus, suggesting that the effects of epinephrine on actinobacteria may be universal.

    Citation: AV Gannesen, MI Schelkunov, RH Ziganshin, MA Ovcharova, MV Sukhacheva, NE Makarova, SV Mart'yanov, NA Loginova, AM Mosolova, EV Diuvenji, ED Nevolina, VK Plakunov. Proteomic and transcriptomic analyses of Cutibacterium acnes biofilms and planktonic cultures in presence of epinephrine[J]. AIMS Microbiology, 2024, 10(2): 363-390. doi: 10.3934/microbiol.2024019

    Related Papers:

    [1] Lin Zhang, Yongbin Ge, Zhi Wang . Positivity-preserving high-order compact difference method for the Keller-Segel chemotaxis model. Mathematical Biosciences and Engineering, 2022, 19(7): 6764-6794. doi: 10.3934/mbe.2022319
    [2] Sunwoo Hwang, Seongwon Lee, Hyung Ju Hwang . Neural network approach to data-driven estimation of chemotactic sensitivity in the Keller-Segel model. Mathematical Biosciences and Engineering, 2021, 18(6): 8524-8534. doi: 10.3934/mbe.2021421
    [3] Thierry Colin, Marie-Christine Durrieu, Julie Joie, Yifeng Lei, Youcef Mammeri, Clair Poignard, Olivier Saut . Modeling of the migration of endothelial cells on bioactive micropatterned polymers. Mathematical Biosciences and Engineering, 2013, 10(4): 997-1015. doi: 10.3934/mbe.2013.10.997
    [4] Shangbing Ai, Zhian Wang . Traveling bands for the Keller-Segel model with population growth. Mathematical Biosciences and Engineering, 2015, 12(4): 717-737. doi: 10.3934/mbe.2015.12.717
    [5] Tong Li, Zhi-An Wang . Traveling wave solutions of a singular Keller-Segel system with logistic source. Mathematical Biosciences and Engineering, 2022, 19(8): 8107-8131. doi: 10.3934/mbe.2022379
    [6] Maghnia Hamou Maamar, Matthias Ehrhardt, Louiza Tabharit . A nonstandard finite difference scheme for a time-fractional model of Zika virus transmission. Mathematical Biosciences and Engineering, 2024, 21(1): 924-962. doi: 10.3934/mbe.2024039
    [7] Qigang Deng, Fugeng Zeng, Dongxiu Wang . Global existence and blow up of solutions for a class of coupled parabolic systems with logarithmic nonlinearity. Mathematical Biosciences and Engineering, 2022, 19(8): 8580-8600. doi: 10.3934/mbe.2022398
    [8] Z. Jackiewicz, B. Zubik-Kowal, B. Basse . Finite-difference and pseudo-spectral methods for the numerical simulations of in vitro human tumor cell population kinetics. Mathematical Biosciences and Engineering, 2009, 6(3): 561-572. doi: 10.3934/mbe.2009.6.561
    [9] Benjamin Wacker, Jan Christian Schlüter . A non-standard finite-difference-method for a non-autonomous epidemiological model: analysis, parameter identification and applications. Mathematical Biosciences and Engineering, 2023, 20(7): 12923-12954. doi: 10.3934/mbe.2023577
    [10] Fadoua El Moustaid, Amina Eladdadi, Lafras Uys . Modeling bacterial attachment to surfaces as an early stage of biofilm development. Mathematical Biosciences and Engineering, 2013, 10(3): 821-842. doi: 10.3934/mbe.2013.10.821
  • Transcriptomic and proteomic analysis were performed on 72 h biofilms of the acneic strain Cutibacterium acnes and planktonic cultures in the presence of epinephrine. Epinephrine predominantly downregulated genes associated with various transporter proteins. No correlation was found between proteomic and transcriptomic profiles. In control samples, the expression of 51 proteins differed between planktonic cultures and biofilms. Addition of 5 nM epinephrine reduced this number, and in the presence of 5 µM epinephrine, the difference in proteomic profiles between planktonic cultures and biofilms disappeared. According to the proteomic profiling, epinephrine itself was more effective in the case of C. acnes biofilms and potentially affected the tricarboxylic acid cycle (as well as alpha-ketoglutarate decarboxylase Kgd), biotin synthesis, cell division, and transport of different compounds in C. acnes cells. These findings are consistent with recent research on Micrococcus luteus, suggesting that the effects of epinephrine on actinobacteria may be universal.



    In this paper, the finite-difference method (FDM) is applied to solve the following two-dimensional (2D) Keller-Segel model [1,2]:

    {ut+(χuv)=dΔu,vt=Δvv+u,(x,y,t)Ω×(0,T]. (1.1)

    The model (1.1) describes the evolutionary process of cell density u(x,y,t) and a chemical stimulus (chemoattractant) concentration v(x,y,t) over a time t and location (x,y), where Ω={(x,y)|ax,yb}R2 is a convex bounded domain and a and b are constants. and Δ are gradient and Laplacian operators, respectively. χ>0 represents the chemotactic sensitivity constant; d>0 denotes the diffusion rate of cells. In addition, the initial conditions related with (1.1) are given as

    u(x,y,0)=u0(x,y),v(x,y,0)=v0(x,y),(x,y)Ω, (1.2)

    and the boundary conditions are assumed to be homogeneous Neumann boundary conditions, that is,

    un=vn=0,(x,y,t)Ω×(0,T], (1.3)

    where Ω is the boundary of Ω and n is the outward normal of Ω. With this condition (1.3), the total mass

    Massu=Ωu(x,y,t)dxdy=Ωu(x,y,0)dxdy

    is conserved as temporal evolution. Besides the above, the model (1.1) has the following form of free energy:

    E(u,v)(t)=Ω(uln(u)+χ2|v|2+χ2v2χuv)dxdy,t>0. (1.4)

    From mathematical analysis, the following equation can be verified by direct calculation of Eq. (1.4),

    E(u,v)t=Ω(u|(ln(u)χv)|2+χ(vt)2)dxdy0,t>0.

    Many early works show that the free energy (1.4) is decreasing over time and is mainly employed to demonstrate the existence of solutions for the chemotaxis system; see [3,4] and the references therein.

    Chemotaxis refers to the directional movement, which includes toward or away from the higher concentrations of cells or microorganisms that are stimulated by chemical substances in the external environment along the gradient directions of the concentration for stimuli. Chemotaxis phenomena play a crucial role in numerous intricate biological evolutions, such as bacterial aggregation, angiogenesis, pattern formation, embryonic development and so on[5]. From as early the 1950 to 1980, the chemotaxis phenomena have been extensively studied by many applied mathematicians and biologists, and a class of models of partial differential systems closely associated with taxis have been proposed; see [1,2,6,7,8,9,10]. Among them the most classical is the above-named system (1.1) (first proposed by Keller and Segel [1,2] in 1970), which simulated the aggregation phenomenon for Amoebae and Dictyostelium, as well as the traveling wave migration phenomenon for Escherichia coli in a capillary filled with nutrients.

    Since the model was proposed, many researchers have systematically analyzed the properties of its solution, including its global existence, asymptotic profile [11], global boundedness [12,13], finite-time blow-up [14,15,16], etc. Particularly, if the initial mass Ωu(x,y)dxdy of the cells in the 2D case satisfies a critical threshold value, its solution will blow up in finite time [7,8,9,14,15,16]. This blow-up denotes a mathematical concept of the bacterial aggregation arising in real biotic environments [7,17,18]. However, it is arduous if we want to obtain the analytical solutions of the model due to the strong nonlinear characteristics of the model itself. Meanwhile, it is still difficult to better numerically capture the blow-up or spike solutions. Therefore, it is desperately needed to establish a high accuracy and more stable numerical method to investigate the properties of the solutions for the chemotaxis system given by (1.1)–(1.3). At the same time, the numerical investigation of the chemotaxis system also facilitates the theoretical exploration of its dynamic behavior.

    In recent years, some numerical methods have been involved via the investigation of the chemotaxis systems [3,4,16,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35]. For example, Saito and Suzuki [19] proposed a conservative numerical scheme by using the FDM to solve a parabolic-elliptic coupling chemotaxis system. Meanwhile, in order to obtain a positivity-preserving scheme under total mass that is conservative, Saito [20] proposed a conservative scheme by using the FDM for the system given by (1.1)–(1.3) with nonlinear diffusion. Xiao et al. [21] derived a semi-implicit scheme by using a characteristic finite element method (FEM) to simulate the blow-up solutions of the chemotaxis system on surfaces, as well as pattern formulation and aggregation phenomenon of bacteria. And, their method has second-order accuracy in L2-norm and H1-norm errors. Epshteyn and Kurganov [22] introduced chemotaxis concentration gradient variables to rewrite the original Keller-Segel model given by (1.1)–(1.3), and they designed an internal penalty discontinuous Galerkin (DG) method for the rewritten system to border on the solution of the original system. Li et al. [23] used a local DG (i.e., LDG) method to modify that in [22], and they obtained the optimal convergence rates based on a specific finite element space before blow-up for the chemotaxis system. They also deduced a positivity-preserving P1 LDG scheme and proved that it is L1-stable. The LDG method was also used to solve the system given by (1.1)–(1.3) in [3], and the energy dissipation with the LDG discretization was proved. Sulman and Nguyen [24] proposed an adaptive moving mesh method by applying an implicit-explicit FEM to solve the system given by (1.1)–(1.3). The method has second-order accuracy in the spatiotemporal directions, and it is noteworthy that the obtained solutions are positive at all time steps if the initial values of the system (1.1) are positive. Qiu et al. [25] proposed a new scheme by using a interface-corrected direct DG method to solve this model (1.1), and their method satisfies the positivity-preserving requiremet without losing third-order accuracy. Based on the gradient flow structure, Shen and Xu [4] proposed a class of numerical schemes for solving the Keller-Segel model given by (1.1)–(1.3). Among them, the first-order accuracy scheme satisfies the mass conservation, bounded positivity, unique solvability and energy dissipation of the original differential equation, while the second-order accuracy scheme satisfies the first three properties. Chen et al. [16] analyzed the error of the numerical scheme proposed in [4], and they deduced the finite-time blow-up of non-radial numerical solutions under certain assumptions. Based on the generalized smoothed particle hydrodynamics meshless method, Dehghan and Abbaszadeh [26] established a second-order-accuracy numerical scheme for some chemotaxis models with the blow-up phenomena during tumor growth, and they numerically simulated the blow-up problems of the original chemotaxis model. Filbet [27] employed the finite volume method (FVM) to approximate the solution of the system given by (1.1)–(1.3), and they simulated the blow-up problems. Chertock and Kurganov [28] developed a center-upward scheme with second-order accuracy by using the FVM based on their previous work in [29] and the blow-up problems of the system (1.1)–(1.3), the pattern formation of bacteria were also simulated. Epshteyn [30] proposed an upwind-difference potentials scheme to solve the problems in complex geometries for the system given by (1.1)–(1.3). In addition, some other numerical methods [31,32,33,34,35], such as the fractional step (or operator splitting) method [31,32], hybrid finite-volume-finite-difference methods [33], generalized FDM [35], etc., were also used to solve the chemotaxis system given by (1.1)–(1.3).

    Although some numerical methods mentioned above for solving the system given by (1.1)–(1.3) can achieve high accuracy in the spatial direction, such as those in [25,30] (third-order accuracy) and [3,33] (fourth-order accuracy), most numerical methods have low accuracy, especially in the temporal direction; see [4,16,19,20,21,22,23,24,26,27,28,29,31,32,34,35]. Meanwhile, many high-accuracy numerical methods are non-compact in space, and the stability conditions are relatively harsh. In other words, if we want to employ these methods to solve real problems, we must take a small time step length to satisfy their stability conditions, which will expend expensive computational time. However, the high-order compact (HOC) difference scheme has attracted many researchers because of its strong advantages, such as fewer computational nodes, small computational errors, better numerical stability and non-complicated boundaries. Meanwhile, the backward differentiation formula with fourth-order accuracy (BDF-4), which is an A-stabilized method and appears first in [36], has been verified to be a feasible method to obtain high accuracy, and it has a relatively large stability condition range [37,38]. In addition, the major difficulty in solving the system comes from the nonlinearities, such as the chemotaxis term (χuv). One issue is that the coupling form will increase the discrete difficulty when we want to obtain the high accuracy and satisfy mass conservative schemes. Another difficulty is positivity preservation since u and v in the system (1.1)-(1.3) have real biological significance in complex chemotaxis phenomenon and cannot have negative values. To acquire the high-accuracy, positivity-preserving numerical solutions for the system given by (1.1)–(1.3), in this work, our purpose was to derive a compact difference scheme to approximate the solutions of the original chemotaxis system, and the scheme has space-time fourth-order accuracy and is stable, as well as positivity-preserving.

    The remainder of this paper is organized as follows. Some preliminary preparations with basic symbols, definitions and theorems are provided in Section 2. In Section 3, we deduce an HOC scheme for the system given by (1.1)–(1.3), and give the computational strategies for the initial time steps and the nonlinear terms. In Section 4, a time advancement algorithm combined with a multigrid method and a positivity-preserving algorithm are proposed. In Section 5, some numerical examples are employed to verify the accuracy, stability, positivity-preserving property, mass conservation and energy dissipation. The finitetime blow-up problems for the chemotaxis system given by (1.1)–(1.3) are simulated by using the proposed method. Finally the conclusion is provided in the end.

    First, the domain {(x,y,t)|ax,yb,0tT} is divided by uniform meshes N2×M, M,NZ+. Denote h=(ba)/N to represent the spatial step size, and τ=T/M stands for the temporal step length. Let xi=a+ih,yj=a+jh and tn=nτ, 0i,jN, 0nM, and mark the mesh points (xi,yj,tn). Figure 1 shows the 2D spatial mesh point stencil.

    Figure 1.  Spatial mesh point stencil.

    Second, we define Ωh={(xi,yj)|0i,jN}; its discrete boundary Ωh={(0,j),(N,j)|0jN}{(i,0),(i,N)|1iN1}; let Ωτ={tn|0nM} and Ωhτ=Ωh×Ωτ. For any mesh function wWhτ={wni,j|0i,jN,0nM} defined on Ωhτ, we have

    wn12i,j=12(wni,j+wn1i,j),δtwn12i,j=1τ(wni,jwn1i,j),δxwni,j=12h(wni+1,jwni1,j),δywni,j=12h(wni,j+1wni,j1),δ2xwni,j=1h2(wni+1,j2wni,j+wni1,j),δ2ywni,j=1h2(wni,j+12wni,j+wni,j1),Δtwni,j=112τ(25wni,j48wn1i,j+36wn2i,j16wn3i,j+3wn4i,j).

    Next, for simplicity, we let pwςp(x,y,t):≜wςp(x,y,t), pZ+, where ς represents x, y or t. And, denote A=[A1,A2,A3,A4,A5]=[25,48,36,16,3], B=[B1,B2,B3]=[1,10,1], S=[S1,S2,S3,S4]=[1,34,13,116]. To obtain a higher accuracy on the boundary described by Eq.(1.3), we employ the Taylor expansion with the Peano remainder to cope with Eq. (1.3). Thus, the following theorem holds.

    Theorem 2.1. Denote Wh={w(xi,yj)|0i,jN}; for mapping w:ΩhWh, we have the following:

    (1) If w(x,y)C5,0([x0,x4]×[y0,yN]), then wx(x0,yj)=112h5k=1Akw(xk1,yj)+O0,j(h4);

    (2) If w(x,y)C5,0([xN4,xN]×[y0,yN]), then wx(xN,yj)=112hN1k=N5ANkw(xk+1,yj)+ON,j(h4);

    (3) If w(x,y)C0,5([x0,xN]×[y0,y4]), then wy(xi,y0)=112h5k=1Akw(xi,yk1)+Oi,0(h4);

    (4) If w(x,y)C0,5([x0,xN]×[yN4,yN]), then wy(xi,yN)=112hN1k=N5ANkw(xi,yk+1)+Oi,N(h4), where the local truncation errors are

    O0,j(h4)=h464k=110k5Skwx5(x0+ksh,yj)(1s)4ds,0jN,ON,j(h4)=h464k=110k5Skwx5(xNksh,yj)(1s)4ds,0jN,Oi,0(h4)=h464k=110k5Skwy5(xi,y0+ksh)(1s)4ds,0iN,Oi,N(h4)=h464k=110k5Skwy5(xi,yNksh)(1s)4ds,0iN.

    Proof. First, we prove (1). We assume that w(x,y)Ck+1,0([xi1,xi+1]×[y0,yN]), and, according to the Taylor expansion with the Peano remainder, we have

    w(xi±h,yj)=kl=0(1)lhll!wxl(xi,yj)+hk+1k!10wxk+1(xi±sh,yj)(1s)kds, (2.1)

    where 1iN1 and 0jN. We suppose that w(x,y)C5,0([x0,x4]×[y0,yN]) and expand w(x1,yj), w(x2,yj), w(x3,yj) and w(x4,yj) at (x0,yj) by using Eq. (2.1) above, that is, we respectively take k=1,2,3,4, and have

    w(xk,yj)=4l=0(kh)ll!wxl(x0,yj)+(kh)52410wx5(x0+ksh,yj)(1s)4ds. (2.2)

    Then, we perform the following the operation: Eq. (2.2)k=1+α× Eq. (2.2)k=2+β× Eq. (2.2)k=3+γ× Eq. (2.2)k=4 for Eq. (2.2), and denote E=[E1,E2,E3,E4]=[1,α,β,γ], where α,β and γ, are constants; then, we can obtain

    4k=1Ekw(xk,yj)=5m=14k=1(kh)m1(m1)!Ekwxm1(x0,yj)+h5244k=1k5Ek10wx5(x0+ksh,yj)(1s)4ds. (2.3)

    Suppose that the coefficients of wx2(x0,yj), wx3(x0,yj) and wx4(x0,yj) in Eq. (2.3) equal to 0; we obtain

    {4k=1k2Ek=0,4k=1k3Ek=0,4k=1k4Ek=0.{α=34,β=13,γ=116.

    Substituting them into Eq. (2.3), denote A=[25,48,36,16,3] and S=[1,34,13,116], and we have

    wx(x0,yj)=112h5k=1Akw(xk1,yj)+O0,j(h4),0jN,

    where O0,j(h4)=h464k=110k5Skwx5(x0+ksh,yj)(1s)4ds. Similarly, (2) and (3) are also easily obtained. The proof is complete.

    According to Theorem 2.1, and by considering wx(x0,yj)=0, wx(xN,yj)=0, wy(xi,y0)=0 and wy(xi,yN)=0, we can easily derive the following fourth-order-accuracy boundary approximation formulas:

    w(x0,yj)125[48w(x1,yj)36w(x2,yj)+16w(x3,yj)3w(x4,yj)],w(xN,yj)125[48w(xN1,yj)36w(xN2,yj)+16w(xN3,yj)3w(xN4,yj)],w(xi,y0)125[48w(xi,y1)36w(xi,y2)+16w(xi,y3)3w(xi,y4)],w(xi,yN)125[48w(xi,yN1)36w(xi,yN2)+16w(xi,yN3)3w(xi,yN4)].

    In addition, we suppose that w(x,y,t)C6,6,5([xi1,xi+1]×[yj1,yj+1]×[tn,tn1]) for 1i,jN1, 1nM, and we denote I=[1,24,81,64] to derive the following truncation errors based on the Taylor expansion with the Peano remainder, that is,

    (Oxx)ni,j(h4)=h424102k=1[13(1s)315(1s)5]wx6(xi+(1)k1sh,yj,tn)ds,(Ot)ni,j(τ4)=τ4610(1μ)44k=1Ikwt5(xi,yj,tnkμτ)dμ,(Ox)ni,j(h4)=h412102k=1[13(1s)314(1s)4]wx5(xi+(1)k1sh,yj,tn)ds,(Ot)n12i,j(τ2)=τ21610(1μ)22k=1wt3(xi,yj,tn12+(1)k1μτ2)dμ,˜On12i,j(τ2)=τ2410(1μ)2k=1wt2(xi,yj,tn12+(1)k1μτ2)dμ.

    Similarly, we can easily derive the expressions of (Oyy)ni,j(h4), (Oy)ni,j(h4), (Ox)n1i,j(h4), (Oy)n1i,j(h4), (Oxx)n12i,j(h4) and (Oyy)n12i,j(h4). To facilitate mathematical analysis below, the following definitions are given.

    Definition 2.2. For any mesh function {wi,j|0i,jN}, the average operators A and B are defined as

    Awi,j={1254k=1Ak+1wk,j,i=0,0jN,1123k=1Bkwi+k2,j,1i,jN1,125N1k=N4ANk+1wk,j,i=N,0jN,Bwi,j={1254k=1Ak+1wi,k,j=0,0iN,1123k=1Bkwi,j+k2,1i,jN1,125N1k=N4ANk+1wi,k,j=N,0iN.

    Definition 2.3. The maximum norm and 2-norm errors are defined as

    ||||=max0i,jN|wMi,jw(xi,yj,tM)|,||||2=h2Ni,j=0[wMi,jw(xi,yj,tM)]2,

    where w(xi,yj,tM) and wMi,j stand for the exact and numerical solutions at the discrete mesh point (xi,yj,tM), respectively.

    Definition 2.4. Denote wmax(t):≜max(x,y)w(x,y,t); the relatively maximum norm and 2-norm errors in the temporal dimension are defined as

    Rel=

    where represents the reference solution.

    Definition 2.5. The convergence rate is defined as

    where stands for or and and are corresponding norm errors which are closely related to and , respectively.

    In this part, an HOC scheme is derived to border on the solution of the Keller-Segel system given by (1.1)–(1.3). To facilitate the numerical analysis, we rewrite an equivalent form for the system given by (1.1)–(1.3), that is,

    where

    Equation (3.1) is a nonlinear advection-diffusion-reaction equation. Then, we define

    on , respectively, and employ the FDM to discretize Eqs. (3.1)–(3.4) for the interior points.

    Now, we focus on Eq. (3.1) at for and ; we have

    (3.5)

    First, the following formulas are applied to discretize the second derivative terms on the right-hand side of Eq. (3.5), i.e.,

    (3.6)

    where , and , and are the central difference operators. Then, for the nonlinear advection terms and on the left-hand side of Eq. (3.5), the following Padé compact schemes [39] are applied to compute them, that is,

    (3.7)

    where and . To be consistent with the accuracy of the interior points above, the following fourth-order boundary schemes [40] are employed, i.e.,

    (3.8)
    (3.9)
    (3.10)
    (3.11)

    where

    Next, we employ the BDF-4 [36] to cope with the temporal derivative on the left-hand side of Eq. (3.5), i.e., for and

    (3.12)

    Substituting Eqs. (3.6)–(3.12) into Eq. (3.5), using the definitions above, combining Eqs. (3.2)–(3.4) and applying Theorem 2.1 for , we have

    where ++, , , and . Then, there exist positive constants , , , and ; it holds that

    Omitting in Eqs. (3.13)–(3.16) and replacing , , and with , , and , respectively, the following HOC scheme for solving the chemotaxis system given by (1.1)–(1.3) for can be obtained:

    Remark 3.1. For the function values containing the unknown time steps in Eqs. (3.17)–(3.20), the following method is used to compute them, where , which stands for the nonlinear cyclic iterative parameter, denotes the approximate convergent solution at the th time step and represents a small amount. That is, for the th time step, the following is performed:

    Approximate using and solve and using Eqs. (3.3), (3.4) and (3.7);

    Compute the chemotaxis terms and using Eqs. (3.7)–(3.11);

    Solve the linear system given by (3.17)–(3.20) for ;

    Set ; if , then proceed to the th time step.

    Remark 3.2. Equations (3.17)–(3.20) demonstrate fully implicit compact scheme including five time steps, and its truncation error is , i.e., it has space-time fourth-order accuracy. It is noteworthy that the compactness here only refers to the spatial direction, since no more than nine mesh points are required for the 2D spatial mesh subdomain. However, it is non-compact in the temporal direction.

    Remark 3.3. For the HOC scheme given by (3.17)–(3.20), besides the values of being known, the values of , and are also needed; then, we can use the scheme given by (3.17)–(3.20) to compute the values of in turn. Thus, in the following part, we discuss the computational strategies for solving , and .

    Next, we consider Eq. (3.1) at for and ; we have

    (3.21)

    First, we employ the Crank-Nicolson (C-N) method to cope with the temporal derivative term in Eq. (3.21), that is,

    Then, we employ the following method for the remaining parts, i.e.,

    where could be and in Eq. (3.21). Then, Eq. (3.21) becomes

    (3.22)

    We use Eq. (3.6) to discretize , , and in the spatial direction of Eq. (3.22), and we use Eqs. (3.7)–(3.11) to compute , , and in the spatial direction. According to Eqs. (3.2)–(3.4) and Theorem 2.1, after rearrangement, the following form can be obtained:

    where . , , , and . Then, there exist positive constants , , , and such that

    Omitting in Eqs. (3.23)–(3.26) and replacing , , and with , , and , respectively, we can get the following scheme for the initial time steps to solve the chemotaxis system given by (1.1)–(1.3) as follows:

    Remark 3.4. The same method in the previous section can be used to compute , , , , and in Eq. (3.27), as described in Remark 3.1.

    Remark 3.5. Equations (3.27)–(3.30) constitute a two-level implicit scheme, and its truncation error is .

    Remark 3.6. Since the C-N method is used to discretize the time derivative term, the scheme given by (3.27)–(3.30) for the initial time steps is unconditionally stable. On the other hand, the formula (3.12) is a -step method for time integration. In this work, only the fourth-order BDF scheme () is used, which is widely employed to solve stiff differential equations [37]. Therefore, according to [36,37,38], the HOC scheme given by (3.17)–(3.20) is -stable.

    As described in Subsections 3.1 and 3.2 above, we use the two-level implicit scheme given by (3.27)–(3.30) to compute , and based on a known . The Richardson extrapolation technique is applied to enhance the accuracy of the scheme given by (3.27)–(3.30) in time and ensure that it is consistent with that of the scheme given by (3.17)–(3.20); the following extrapolation formula is used, that is,

    (3.31)

    Here, stands for the extrapolated solution at the th time step. represents the computed solution obtained by using the scheme given by (3.27)–(3.30) with the temporal step length , whereas denotes the computed solution obtained by using the scheme given by (3.27)–(3.30) with .

    Remark 3.7. Equation (3.31) can extrapolate the accuracy in the temporal direction from second to fourth order.

    In this part, based on the compact difference schemes proposed above, we will establish the corresponding numerical algorithm. On one hand, because the slow convergence speed of the classical iterative method (for instance, Gauss–Seidel, Jacobi or successive over-relaxation iterations) used to solve the algebraic systems arise out of the fully implicit scheme at each time step, we aim to employ a multigrid method to accelerate the convergence speed. On the other hand, we want to establish a positivity-preserving algorithm to obtain the non-negative solutions of and without losing the fourth-order accuracy. Finally, we structure a time advancement algorithm based on these algorithms above to solve the system given by (1.1)–(1.3).

    Since the idea of the multigrid method was proposed in the 1930s, until Professor Brandt published his pioneering work [41] in 1977, this method was increasingly applied to solve engineering and technical problems. Up to now, the multigrid method has been widely applied in numerous disciplines and engineering fields, such as computational fluid dynamics and computational biology. The multigrid method has been theoretically proved to be a first-rank numerical computational method for linear elliptic problems [41,42]. Its convergence rate is independent of the mesh scale, and the computational speed of the existing computational program using the classical iterative method can be increased by employing the multigrid method with to orders of magnitude, which is especially suitable for application in large-scale engineering numerical computational problems. Nowadays, the multigrid method is also used more and more to solve non-elliptic problems to speed up the convergence of each time step and satisfy the needs of practical problems [43,44].

    The multigrid method is implemented by using a cycling algorithm. A multigrid cycle includes three elements: relaxation, restriction and interpolation operators. We use alternating direction line Gauss–Seidel (ALGS) [45] iteration for the relaxation operator, as well as half-weighted and fully weighted restriction operators and a bilinear interpolation operator [42]. As the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are both nonlinear, we employ a multigrid full approximation scheme (FAS) [42,44] to accommodate nonlinearities. For simplicity, we formally denote the algebraic equations arising from Eqs. (3.27)–(3.30) and (3.17)–(3.20) at each time step as

    (4.1)

    The multigrid FAS algorithm is implemented recursively. Here, we only give the following two-level FAS cycle algorithm based on Eq. (4.1):

    Algorithm 4.1: Two-level FAS cycle algorithm
    ;

    Remark 4.1. Here, and are the difference operators; (half-weighted) and (fully weighted) are the restriction operators, which are used to restrict the approximate value and residual from the fine mesh level to coarse mesh level ; and, is the interpolation operator, which is used to transfer the error correction from the coarse to fine mesh levels. is the residual at the fine mesh level, and is the residual at the adjacent coarse mesh level; and represent the pre-smooth and post-smooth numbers, respectively.

    Remark 4.2. In this work, we employ fully multi-level algorithms for all calculations. For instance, if the finest grid is , an eight-level algorithm is used. Under such conditions, the coarsest grid level only has grids.

    When the schemes given by (3.27)–(3.30) and (3.17)–(3.20) are employed to approximate the solution of the chemotaxis system, it may generate negative values where blow-up occurs. To obtain the non-negative values of and for all time levels, we define the following average solution for each time step:

    (4.2)

    where , , and . Similar to the numerical integration formula of the one-variable function (i.e., Simpson formula), we structure the numerical integration of the two-variable function on for each time step, i.e., for ,

    (4.3)

    Substituting Eq. (4.3) into Eq. (4.2) at each time step, we obtain the following average solutions , that is, for and , we have

    (4.4)

    Next, a positivity-preserving limiter (PPL) [25,46] is employed to eliminate the negative values of , , , that is,

    (4.5)

    where For and , we can get the positivity-preserving solution without losing the fourth-order accuracy obtained by the proposed scheme. Finally, we establish the following algorithm for the th time step:

    Algorithm 4.2: Positivity-preserving algorithm
    --

    On the basis of the derivation process described in Section 3 and the proposed algorithms in Subsections 4.1 and 4.2, we employ the difference schemes given by (3.27)–(3.30) and (3.17)–(3.20) to approximate the solution of the original system given by (1.1)–(1.3), which mainly includes the following three steps: first, provide the initial value of ; second, we employ Eqs. (3.7)–(3.11) and (3.27)–(3.31) to compute the values of , and ; finally, we employ Eqs. (3.7)–(3.11) and (3.17)–(3.20) to compute the value of . As (3.27)–(3.30) and (3.17)–(3.20) involve nonlinear parts, we establish the following algorithm with a nonlinear iteration tactic to solve them.

    Algorithm 4.3: Time advancement algorithm
    ;
      
        ;
        ;
        ;
        ;
        
        
        
        ;
        ;
        ;
        ;
        
        
          ;
        
      
        
      
        
        
        ;
        ;
        ;
        
        ;
        
       
      
       

    In this part, we give several numerical experiments to test the various properties of the proposed scheme and algorithm when solving the system given by (1.1)–(1.3). First, in Subsection 5.1 below, the accuracy and stability of the proposed scheme and algorithm in the absence of PPL and in the presence of a PPL when solving the chemotaxis system are tested. Second, in Subsection 5.2 below, we simulate the blow-up phenomena for the system given by (1.1)–(1.3), compare the obtained results from the proposed scheme and algorithm before and after using the positivity-preserving algorithm and verify the mass conservation and energy dissipation of the proposed method.

    We consider the following general type of Eq. (1.1) with source terms to test the accuracy, that is,

    We use the periodic boundary conditions, and its analytical solution is By observing the expression of this analytical solution, we can see that the global solutions of this problem are positive values in the entire physical domain. Therefore, in this computation, we will use the difference scheme proposed in this work to approximate the solution for this problem in the absence of a PPL.

    In Table 1, we take and test the convergence of the HOC scheme at the final time for Problem 5.1.1. From this table, the fourth-order accuracy is obtained by using the proposed scheme. The obtained results for the LDG method in [3] (the time step size is used to solve this problem) represent third-order accuracy, which shows the advantage of our method, that is, we can obtain higher accuracy with a larger time step. Table 2 displays the convergence of the HOC scheme at and for Problem 5.1.1, where . Table 3 shows the errors of and for different step ratios for Problem 5.1.1, where and . From these tables, we obtain that the proposed method still converges at a fourth-order rate when computing long times and large step ratios, which indicates that our method has better stability.

    Table 1.  Convergence of the HOC scheme for Problem 5.1.1 at .
    LDG scheme [3] HOC scheme
    Rate Rate Rate Rate
    8 5.79E-03 2.53E-02 9.415E-03 3.974E-02
    16 8.79E-04 2.69 3.87E-03 2.71 3.359E-04 4.81 1.379E-03 4.85
    32 1.19E-04 2.91 5.18E-04 2.90 1.912E-05 4.13 7.711E-05 4.16
    64 1.49E-05 3.00 6.57E-05 2.98 1.112E-06 4.10 4.431E-06 4.12

     | Show Table
    DownLoad: CSV
    Table 2.  Convergence of the HOC scheme for Problem 5.1.1 with .
    Rate Rate Rate Rate
    Cell density
    20 1.472E-03 6.477E-03 5.428E-06 2.415E-05
    40 1.046E-04 3.81 4.621E-04 3.81 2.289E-07 4.57 1.018E-06 4.57
    80 6.404E-06 4.03 2.832E-05 4.03 1.840E-08 3.64 5.259E-08 4.27
    160 3.902E-07 4.04 1.726E-06 4.04 6.967E-10 4.72 3.043E-09 4.11
    Chemoattractant concentration
    20 1.459E-03 6.488E-03 5.428E-06 2.415E-05
    40 1.041E-04 3.81 4.629E-04 3.81 2.289E-07 4.57 1.018E-06 4.57
    80 6.382E-06 4.03 2.837E-05 4.03 1.840E-08 3.64 5.259E-08 4.27
    160 3.889E-07 4.04 1.729E-06 4.04 6.967E-10 4.72 3.043E-09 4.11

     | Show Table
    DownLoad: CSV
    Table 3.  errors of the HOC scheme for Problem 5.1.1 with .
    0.4 2.8366512125E-05 2.8485324431E-05
    0.8 4.2917900191E-04 4.2928242249E-04
    1.6 5.7178346316E-03 5.7179298620E-03
    0.8 2.8033654005E-05 2.8041090035E-05
    1.6 4.1876675539E-04 4.1877318296E-04
    3.2 5.5523910128E-03 5.5523969211E-03
    1.6 2.7710430090E-05 2.7710912582E-05
    3.2 4.1352388720E-04 4.1352430255E-04
    6.4 5.4716502752E-03 5.4716506553E-03

     | Show Table
    DownLoad: CSV

    Next, we consider the following general type of Eq. (1.1) with two additional fluxes[24]:

    where , and the homogeneous Neumann boundary condition (1.3) is applied.

    Case 1: The additional fluxes are , , with the initial condition , and the analytical solution[23,24]

    Case 2: The additional fluxes are , , with the initial condition , and the analytical solution .

    Table 4 shows the convergence of the HOC scheme when we take the time step size at when using the proposed method to solve Problem 5.1.2 Case 1: the obtained results are compared with those in [24]. According to Table 4, the results obtained via the proposed method in the absence of a PPL can converge at the fourth-order rate. After using the positivity-preserving algorithm, although the errors obtained via the proposed scheme in the presence of a PPL decreases slightly, it still converges at the fourth-order rate and is better than the results in [24], as it is only second-order accuracy in [24]. In Tables 5 and 6, the results computed by using the proposed method for Problem 5.1.2 are given respectively. And, the convergence of the HOC scheme in the presence and absence of a PPL are compared when , and . By observing these tables, it can be found that the computed results for the proposed method can achieve the fourth-order accuracy in both cases, including with and without a PPL. Finally, we take the spatial meshes and , respectively, at the final time . By using the proposed method to solve Problem 5.1.2 Case 1, we obtain the norm errors computed with and without a PPL for different step ratios . The computed results are listed in Table 7. On the basis of this table, with the increase of the step ratio , the proposed method can converge well regardless of PPL existence, which also reflects that the proposed HOC difference scheme has better stability.

    Table 4.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 with .
    Ref. [24] HOC scheme
    Without PPL With PPL
    N Rate Rate Rate
    10 1.02E-01 9.275E-03 2.170E-01
    20 2.74E-02 1.90 6.000E-04 3.95 2.870E-02 2.92
    40 6.30E-03 2.12 2.624E-05 4.52 2.469E-03 3.54
    80 1.45E-03 2.12 1.593E-06 4.04 1.476E-04 4.06

     | Show Table
    DownLoad: CSV
    Table 5.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 when .
    Without PPL With PPL
    Rate Rate Rate Rate
    Cell density
    16 4.670E-03 1.022E-02 4.846E-03 1.424E-02
    32 1.009E-04 5.53 2.447E-04 5.38 2.684E-04 4.17 7.271E-04 4.29
    64 1.686E-06 5.90 5.902E-06 5.37 2.048E-05 3.71 4.974E-05 3.87
    128 7.683E-08 4.46 2.329E-07 4.66 1.508E-06 3.76 3.539E-06 3.81
    Chemoattractant concentration
    16 3.154E-03 8.025E-03 3.457E-03 1.217E-02
    32 9.167E-05 5.10 2.391E-04 5.07 2.769E-04 3.64 7.280E-04 4.06
    64 2.233E-06 5.36 6.687E-06 5.16 2.129E-05 3.70 5.074E-05 3.84
    128 4.873E-08 5.52 1.986E-07 5.07 1.558E-06 3.77 3.574E-06 3.83

     | Show Table
    DownLoad: CSV
    Table 6.  Convergence of the HOC scheme for Problem 5.1.2 Case 1 when .
    Without PPL With PPL
    Rate Rate Rate Rate
    16 2.693E-03 1.418E-02 4.424E-03 2.430E-02
    32 6.517E-05 5.37 3.530E-04 5.33 2.106E-04 4.39 1.172E-03 4.37
    64 1.631E-06 5.32 9.106E-06 5.28 1.331E-05 3.98 7.517E-05 3.96
    128 4.703E-08 5.12 2.672E-07 5.09 8.955E-07 3.89 5.098E-06 3.88

     | Show Table
    DownLoad: CSV
    Table 7.  errors of the HOC scheme for Problem 5.1.2 Case 1 with .
    Without PPL With PPL
    0.4 6.4023871303E-06 7.2616742789E-06 6.1464555035E-05 6.2581655559E-05
    0.8 1.2105436765E-05 1.1596839477E-05 3.5050251455E-05 3.5272598246E-05
    1.6 1.4457697771E-04 1.4260648214E-04 1.4620731799E-04 1.4406453705E-04
    0.8 9.4507876168E-07 8.4039918076E-07 2.4504429125E-06 2.4129793104E-06
    1.6 1.4859666027E-05 1.4234907841E-05 1.4967935575E-05 1.4314292650E-05
    3.2 1.5091392012E-04 1.4863144808E-04 1.5093408941E-04 1.4864245029E-04
    1.6 1.0545853150E-06 9.9391781457E-07 1.0618094946E-06 9.9866861508E-07
    3.2 1.5143113254E-05 1.4499368302E-05 1.5145042129E-05 1.4500157168E-05
    6.4 1.5233575478E-04 1.4973374665E-04 1.5233626588E-04 1.4973395499E-04

     | Show Table
    DownLoad: CSV

    Table 8 lists the numerical convergence of HOC scheme Problem 5.1.2 Case 2 in the absence and presence of a PPL, where and , and it compares the computed results with those obtained via the method in [25]. The results in [25] have third-order accuracy before and after using a PPL, while the results computed via the method in this work converge at the fourth-order rate for both cases, including without and with a PPL, and the , errors are better than those in [25]. This also reflects the superiority of the HOC difference scheme proposed in this work.

    Table 8.  Convergence of the HOC scheme for Problem 5.1.2 Case 2 at .
    Ref. [25] HOC scheme
    Rate Rate Rate Rate
    Without PPL
    20 3.65E-03 3.04E-03 6.189E-04 9.217E-04
    40 4.55E-04 3.00 3.67E-04 3.05 3.554E-05 4.12 7.089E-05 3.70
    80 5.68E-05 3.00 4.58E-05 3.00 2.535E-06 3.81 4.677E-06 3.92
    160 7.02E-06 3.01 5.77E-06 2.99 1.741E-07 3.86 3.077E-07 3.93
    With PPL
    20 3.55E-02 3.62E-02 6.337E-04 1.019E-03
    40 4.37E-03 3.02 4.54E-03 3.00 3.788E-05 4.06 7.260E-05 3.81
    80 5.40E-04 3.02 5.61E-04 3.02 2.562E-06 3.89 4.706E-06 3.95
    160 6.74E-05 3.01 6.92E-05 3.02 1.745E-07 3.88 3.082E-07 3.93

     | Show Table
    DownLoad: CSV

    In the 2D case of the system given by (1.1)–(1.3), when the initial mass satisfies certain critical values, the solutions will blow up in a finite time. To capture the blow-up time of the solutions, referring to [47], the following adaptive technology is used to obtain the optimal time step sizes, i.e.,

    (5.1)

    The computation is terminated if , where denotes the maximum order of magnitude of cell density when the actual problem experiences a blow up. And, is regarded as the approximation of the blow-up time .

    In a rectangular region , we first take and test the initial-boundary value problem (IBVP) [23,24,25] for the chemotaxis system given by (1.1)–(1.3) with the homogeneous Neumann boundary conditions given by (1.3), where its initial conditions are as follow:

    The solution of this problem will blow up in finite time at the center of this rectangular region.

    (1) Finite-time blow up. Figure 2 plots the numerical solutions of cell density as obtained by using the proposed method in the presence of a PPL, to solve the IBVP 5.2.1 at , , and , where the initial time step size and the fixed spatial mesh numbers are . In [22], the blow-up time is approximately equal to , which is verified by the authors. Based on observation, we can expect to see that Problem 5.2.1 displays blow-up phenomena around at the center of the rectangular region . It shows that the maximum peak values of the cell density gradually increase over time and present a very sharp aiguille structure at the central zone. We find that is strictly positive during time evolution by using the proposed method in the presence of a PPL. Figure 3 displays the projections of on the plane for Problem 5.2.1 at two pre-blow-up times and , with and . They intuitively show that the proposed method in this work can achieve positivity-preserving capability.

    Figure 2.  Cell density with PPL for Problem 5.2.1.
    Figure 3.  Projection for on the plane for Problem 5.2.1.

    (2) Non-negativity. To highlight the performance of positivity-preserving capability for our method in the presence of a PPL, we employ the proposed scheme in the absence and presence of a PPL to solve the IBVP 5.2.1 in two different space grids and at pre-blow-up times and , where . In Figure 4, we can see that has negative values in the absence of a PPL, and that all negative values are eliminated after the PPL is used. At the same time, the negative values decrease with the increase of spatial grid numbers, and the solutions become smoother. However, for the same number of grids, the oscillation of the blow-up area increases with time.

    Figure 4.  One-dimensional profile of along for Problem 5.2.1.

    (3) Influence of chemotaxis sensitivity coefficient . We employ the proposed method to compute the blow-up solution in a finite time for this IBVP 5.2.1 under different chemotaxis sensitivity coefficients and at the same pre-blow-up time . The results are shown in Figure 5, where and . We find that, for different values of , the cell density achieves negative values in the absence of a PPL. After using the PPL, the algorithm effectively eliminates the negative values. At the same time, the maximum peak value when is one order of magnitude higher than that when . It can be seen that the peak value is gradually increased over . Therefore, the sensitivity intensity also affects the value of and the blow-up degree of Problem 5.2.1.

    Figure 5.  One-dimensional profile of along for Problem 5.2.1 with .

    (4) Chemoattractant concentration. In Figure 6, we plot the numerical solutions of chemoattractant concentration for Problem 5.2.1 when , and , and we have compared the computed results in the absence and presence of a PPL. Observing Figures 4b and 6, we find that, although we applied the same parameters in the absence of a PPL, achieves negative values, but does not have negative values, and that the numerical solutions in the presence of a PPL is highly consistent with those in the absence of a PPL. In summary, our method can effectively eliminate these negative values in this computation. It is consistent with the results obtained in the absence of a PPL if the obtained solutions do not have negative values. Meanwhile, Problem 5.2.1 always experiences a finite-time blow-up at the center of the rectangular region under this initial condition.

    Figure 6.  Chemoattractant concentration for Problem 5.2.1 with .

    (5) Accuracy and mass conservation. Since the analytical solution for this IBVP 5.2.1 is very difficult to obtain, to verify the accuracy and mass conservation, we take the fine mesh as the reference solution. The convergence of the HOC scheme in the presence of a PPL are tested by using the relative errors between the maximum value of and the reference solutions ; the results are shown in Table 9. We find that can converge at a fourth-order relative error convergence rate. Meanwhile, we plot , which is the time evolution of the maximum value of , on the left of Figure 7, and the numerical energy with time on the right of Figure 7, using the proposed scheme in the presence of a PPL to solve this IBVP 5.2.1 with PPL at , , . We find that the energy is strictly positive and monotonically decreasing with the evolution of time, which verifies the energy dissipation nature of the original problem. In addition, in Table 10, we give the discrete masses of and as obtained by using the proposed scheme in the absence and presence of a PPL to solve this IBVP 5.2.1 at different times with . According to this table, we can find that the proposed scheme in the absence of a PPL well verifies the mass conservation. After using the PPL, it is still conserved before blow up, and the mass conservation is slightly affected near the blow-up time.

    Table 9.  Relative , errors and convergence rates of for Problem 5.2.1 with .
    N Rate Rate
    100 3.204E-02 2.104E-02
    200 1.864E-03 4.10 1.197E-03 4.14
    300 3.290E-04 4.28 2.134E-04 4.25
    400 7.872E-05 4.97 5.902E-05 4.47
    500 3.332E-05 3.85 2.154E-05 4.52

     | Show Table
    DownLoad: CSV
    Table 10.  and of the HOC scheme for Problem 5.2.1 at time with .
    Without PPL With PPL
    31.4159265323 31.4156968041 31.4159265323 31.4156968041
    31.4159265336 31.4157188735 31.4159265336 31.4157188735
    31.4159265339 31.4157222761 31.4159265339 31.4157222761
    31.4159265341 31.4157255274 31.4159265341 31.4157255274
    31.4159265344 31.4157287900 31.4159265344 31.4157287900
    31.4159265348 31.4157319661 31.4159265348 31.4157319661
    31.4159265354 31.4157381680 31.4292462953 31.4157381680
    31.4159265366 31.4157472806 37.7315206418 31.4157952866
    31.4159265376 31.4157533114 48.2548641463 31.4159976499

     | Show Table
    DownLoad: CSV
    Figure 7.  Left: Plots of : the time evolution of the maximum value for ; Right: Plots of the free energy with time for Problem 5.2.1 with the PPL at and .

    In the section, we consider the following IBVP [24] for the chemotaxis system given by (1.1)–(1.3) with the boundary conditions given by (1.3) for the rectangular region . Its initial condition is given as follows:

    In the computation, we take and employ the proposed scheme in the presence of a PPL to simulate the temporal evolution of and . The solution for this Problem 5.2.2 will blow up in finite time [8,9,24] at the corner of the rectangular region . Next, we take the initial time step size and fix a mesh of to compute them.

    Figures 8 and 9 respectively display the evolutionary processes of and over time . It can be seen in these figures that the maximum value of gradually moves toward the corner of the rectangular region in the form of a Gaussian profile. When it is close to the corner of the rectangular region, the value of increases rapidly, finally blowing up in a finite time. At the same time, also increases rapidly at the corner of the rectangular region.

    Figure 8.  Evolution of over time for Problem 5.2.2.
    Figure 9.  Evolution of over time for Problem 5.2.2.

    In this work, first, a fourth-order compact difference scheme for solving a 2D Keller-Segel model was derived and the computing strategies for the initial time steps and the nonlinear chemotaxis terms has been presented. Then, a multigrid algorithm was established to improve the computational efficiency and the binary function integration method has been used to establish the positivity-preserving algorithm. Third, the accuracy and reliability of the proposed method in the absence and presence of a PPL have been respectively verified by several numerical examples with flux source terms. And, we simulated the blow-up phenomena for the chemotaxis system in the center and corner of the rectangular region, respectively. Finally, the energy dissipation and mass conservation were also verified by numerical experiments. In summary, compared with the classical high-order numerical methods for solving the chemotaxis system in the literature, this proposed method has the following advantages:

    ● The proposed HOC scheme is compact in the spatial directions and fully implicit, as no more than nine mesh points are required for the 2D spatial mesh subdomain.

    ● The truncation error of the proposed HOC scheme is , i.e., it has space-time fourth-order accuracy. Thus, we can take large time step sizes to approximate the solution of this chemotaxis system.

    ● The proposed numerical algorithm can effectively filter the negative values in the solution process to ensure that the values of cell density are not negative at all time steps, while maintaining the original fourth-order accuracy without loss.

    ● The computed results show that the proposed method is more accurate than most such schemes reported in the literature.

    Nevertheless, our method also has some disadvantages. First, the new scheme has five layers in time. When it is used to compute the actual problem, three start-up time steps are required, which is slightly more troublesome for programming. Second, the proposed positivity-preserving algorithm would slightly affect its mass conservation when the actual problem blows up. The above deficiencies are also the driving force for our next work. At the same time, it is also possible to extend our approach to solve more complex chemotaxis and haptotaxis systems, which will be incorporated into our next work.

    Support of the study was received from the National Natural Science Foundation of China (12161067, 12001015, 12261067), National Natural Science Foundation of Ningxia, China (2022AAC02023) and National Youth Top-Notch Talent Support Program of Ningxia, China.

    All authors assert no conflict of interest.

    The data that support the findings of this study are available from the authors upon reasonable request.


    Acknowledgments



    Transcriptome sequencing and analysis were performed in the Genomics Core Facility of the Skolkovo Institute of Science and Technology (https://www.skoltech.ru/research/en/gcf-2/ (accessed on March1, 2024). qRT–PCR was performed using the scientific equipment of the Core Research Facility ‘Bioengineering’ (Research Center of Biotechnology RAS). A.V. Gannesen and coauthors sincerely thank V.V. Sorokin and A.V. Mulyukin (Core Facility “UNIQEM collection” of Research Center of Biotechnology RAS) for their kind assistance and for providing liquid nitrogen for this study. The authors are grateful to the staff of the Laboratory of Molecular Ecology and Phylogenomics of Bacteria of the Research Center of Biotechnology RAS (the head Svetlana N. Dedysh) and, personally to Svetlana E. Belova, Igor. Yu. Oshkin, Kirill K. Miroshnikov, and Anastasiya A. Ivanova for their kind assistance and giving access to the FastPrep disintegrator. Authors cordially thank Maya M. Polovitskaya for her kind assistance in the editing of the manuscript.

    Funding



    The work is supported by the RSCF, the grant nº 19-74-10071. The work of V.K. Plakunov is supported by the Ministry of Science and Higher Education of Russian Federation.

    Conflict of interests



    Authors declare no conflict of interests.

    Data availability statement



    The transcriptomic data are deposited at the Sequence Read Archive (SRA) database under the identifiers SRR27884912 and SRR27884920. The mass spectrometry proteomics data have been deposited at the ProteomeXchange Consortium via the PRIDE partner repository [57][58] with the data set identifier PXD049410.

    Author contributions



    Conceptualization, A.V.G.; investigation, M.A.O., M.I.S., R.H.Z., N.A.L., A.M.M., N.E.M., S.V.M., E.D.N., E.D.V. and A.V.G.; validation, M.A.O., M.I.S., S.V.M., R.H.Z., N.A.L., A.M.M., N.E.M., S.V.M., E.D.N., V.K.P., and A.V.G.; formal analysis, M.A.O., M.I.S., R.H.Z., and A.V.G.; writing—original draft, M.A.O., M.I.S., S.V.M., R.H.Z., E.D.N., N.E.M., N.A.L., A.M.M., E.D.V., and A.V.G.; data curation, M.A.O., M.I.S., R.H.Z., N.E.M., and A.V.G.; writing—review and editing, M.A.O., M.I.S., R.H.Z., S.V.M., V.K.P. and A.V.G.; visualization, M.A.O., M.I.S., R.H.Z., S.V.M., E.D.N., N.A.L., A.M.M., and A.V.G.; supervision, V.K.P. and A.V.G.; resources, A.V.G.; funding acquisition, A.V.G.; project administration, A.V.G. All authors have read and agreed to the published version of the manuscript.

    [1] Lyte M (2014) Microbial endocrinology and the microbiota-gut-brain axis. Microbial Endocrinology: The MicrobiotaGut-Brain Axis in Health and Disease. New York: Springer 3-24. https://doi.org/10.1007/978-1-4939-0897-4_1
    [2] Luqman A (2023) The orchestra of human bacteriome by hormones. Microb Pathog 106125. https://doi.org/10.1016/j.micpath.2023.106125
    [3] Niu L, Gao M, Wen S, et al. (2023) Effects of catecholamine stress hormones norepinephrine and epinephrine on growth, antimicrobial susceptibility, biofilm formation, and gene expressions of enterotoxigenic Escherichia coli. Int J Mol Sci 24: 15646. https://doi.org/10.3390/ijms242115646
    [4] Reading NC, Rasko D, Torres AG, et al. (2009) A transcriptome study of the QseEF two-component system and the QseG membrane protein in enterohaemorrhagic Escherichia coli O157: H7. Microbiology 156: 1167. https://doi.org.10.1099/mic.0.033027-0
    [5] Moreira CG, Sperandio V (2010) The epinephrine/norepinephrine/autoinducer-3 interkingdom signaling system in Escherichia coli O157: H7. Microbial Endocrinology: Interkingdom Signaling in Infectious Disease and Health. New York: Springer 213-227. https://doi.org/10.1007/978-1-4419-5576-0_12
    [6] Lopes JG, Sourjik V (2018) Chemotaxis of Escherichia coli to major hormones and polyamines present in human gut. ISME J 12: 2736-2747. https://doi.org/10.1038/s41396-018-0227-5
    [7] Qin T, Xie J, Xi B, et al. (2021) Integrative transcriptomic and proteomic analyses of pathogenic Aeromonas hydrophila in response to stress hormone norepinephrine. Aquac Res 53: 1693-1705. https://doi.org/10.1111/are.15700
    [8] Cambronel M, Nilly F, Mesguida O, et al. (2020) Influence of catecholamines (epinephrine/norepinephrine) on biofilm formation and adhesion in pathogenic and probiotic strains of Enterococcus faecalis. Front Microbiol 11: 1501. https://doi.org/10.3389/fmicb.2020.01501
    [9] Cambronel M, Tortuel D, Biaggini K, et al. (2019) Epinephrine affects motility, and increases adhesion, biofilm and virulence of Pseudomonas aeruginosa H103. Sci Rep 9: 20203. https://doi.org/10.1038/s41598-019-56666-7
    [10] Intarak N, Muangsombut V, Vattanaviboon P, et al. (2014) Growth, motility and resistance to oxidative stress of the melioidosis pathogen Burkholderia pseudomallei are enhanced by epinephrine. Pathog Dis 72: 24-31. https://doi.org/10.1111/2049-632X.12181
    [11] Perraud Q, Kuhn L, Fritsch S, et al. (2022) Opportunistic use of catecholamine neurotransmitters as siderophores to access iron by Pseudomonas aeruginosa. Environ Microbiol 24: 878-893. https://doi.org/10.1111/1462-2920.15372
    [12] Beasley FC, Marolda CL, Cheung J, et al. (2011) Staphylococcus aureus transporters Hts, Sir, and Sst capture iron liberated from human transferrin by Staphyloferrin A, Staphyloferrin B, and catecholamine stress hormones, respectively, and contribute to virulence. BMC Microbiol 11: 199. https://doi.org/10.1128/IAI.00117-11
    [13] Boukerb AM, Cambronel M, Rodrigues S, et al. (2021) Inter-kingdom signaling of stress hormones: sensing, transport and modulation of bacterial physiology. Front Microbiol 12: 690942. https://doi.org/10.3389/fmicb.2021.690942
    [14] Axelrod J, Gordon Ayhitby L, Hertting G, et al. (1961) Studies on the metabolism of catecholamines. Circ Res 9: 715-719. https://doi.org/10.1161/01.RES.9.3.715
    [15] Notelovitz M, Funk S, Nanavati N, et al. (2002) Estradiol absorption from vaginal tablets in postmenopausal women. Obstet Gynecol 99: 556-562. https://doi.org/10.1016/S0029-7844(01)01385-0
    [16] Leis S, Drenkhahn S, Schick C, et al. (2004) Catecholamine release in human skin—a microdialysis study. Exp Neurol 188: 86-93. https://doi.org/10.1016/j.expneurol.2004.03.013
    [17] Stradford AF, Zoutman DE, Dick E, et al. (2002) Effect of lidocaine and epinephrine on Staphylococcus aureus in a guinea pig model of surgical wound infection. Plast Reconstr Surg 110: 1275-1279. https://doi.org/10.1097/01.PRS.0000025427.86301.8A
    [18] Lyte M, Freestone PPE, Neal CP, et al. (2011) Stimulation of Staphylococcus epidermidis growth and biofilm formation by catecholamine inotropes. Lancet 361: 130-135. https://doi.org/10.1016/S0140-6736(03)12231-3
    [19] Mart'yanov SV, Botchkova EA, Plakunov VK, et al. (2021) The impact of norepinephrine on mono-species and dual-species staphylococcal biofilms. Microorganisms 9: 820. https://doi.org/10.3390/microorganisms9040820
    [20] Gannesen AV, Schelkunov MI, Geras'kina OV, et al. (2021) Epinephrine affects gene expression levels and has a complex effect on biofilm formation in Micrococcus luteus strain C01 isolated from human skin. Biofilm 3: 100058. https://doi.org/10.1016/j.bioflm.2021.100058
    [21] Gannesen AV, Ziganshin RH, Zdorovenko EL, et al. (2022) Epinephrine extensively changes the biofilm matrix composition in Micrococcus luteus C01 isolated from human skin. Front Microbiol 13: 1003942. https://doi.org/10.3389/fmicb.2022.1003942
    [22] Gannesen AV, Ziganshin RH, Ovcharova MA, et al. (2023) Epinephrine affects ribosomes, cell division, and catabolic processes in Micrococcus luteus skin strain C01: revelation of the conditionally extensive hormone effect using orbitrap mass spectrometry and proteomic analysis. Microorganisms 11: 2181. https://doi.org/10.3389/fmicb.2022.1003942
    [23] Fournière M, Latire T, Souak D, et al. (2020) Staphylococcus epidermidis and Cutibacterium acnes: two major sentinels of skin microbiota and the influence of cosmetics. Microorganisms 8: 1752. https://doi.org/10.3390/microorganisms8111752
    [24] Akaza N, Takasaki K, Nishiyama E, et al. (2022) The microbiome in comedonal contents of inflammatory acne vulgaris is composed of an overgrowth of Cutibacterium Spp. and other cutaneous microorganisms. Clin Cosmet Investig Dermatol 21: 2003-2012. https://doi.org/10.2147/CCID.S379609
    [25] Jahns AC, Alexeyev OA (2014) Three dimensional distribution of Propionibacterium acnes biofilms in human skin. Exp Dermatol 23: 687-689. https://doi.org/10.1111/exd.12482
    [26] Chiu A, Chon SY, Kimball AB (2003) The response of skin disease to stress: changes in the severity of acne vulgaris as affected by examination stress. Arch Dermatol 139: 897-900. https://doi.org/10.1001/archderm.139.7.897
    [27] Yosipovitch G, Tang M, Dawn AG, et al. (2007) Study of psychological stress, sebum production and acne vulgaris in adolescents. Acta Derm Venereol 87: 135-139. https://doi.org/10.2340/00015555-0231
    [28] Jusuf NK, Putra IB, Sutrisno AR (2021) Correlation between stress scale and serum substance P level in acne vulgaris. Int J Gen Med 1: 681-686. https://doi.org/10.2147/IJGM.S294509
    [29] Boyanova L (2017) Stress hormone epinephrine (adrenaline) and norepinephrine (noradrenaline) effects on the anaerobic bacteria. Anaerobe 44: 13-19. https://doi.org/10.1016/j.anaerobe.2017.01.003
    [30] Borrel V, Thomas P, Catovic C, et al. (2019) Acne and stress: Impact of catecholamines on Cutibacterium acnes. Front Med 6: 155. https://doi.org/10.3389/fmed.2019.00155
    [31] Ovcharova MA, Schelkunov MI, Geras'kina OV, et al. (2023) C-Type natriuretic peptide acts as a microorganism-activated regulator of the skin commensals Staphylococcus epidermidis and Cutibacterium acnes in dual-species biofilms. Biology 12: 436. https://doi.org/10.3390/biology12030436
    [32] Illumina: bcl2fastq and bcl2fastq2 Conversion Software (2024). Available from: https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html
    [33] Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30: 2114-2120. https://doi.org/10.1093/bioinformatics/btu170
    [34] Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25: 1754-1760. https://doi.org/10.1093/bioinformatics/btp324
    [35] Patro R, Duggal G, Love MI, et al. (2017) Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods 14: 417-419. https://doi.org/10.1038/nmeth.4197
    [36] Camacho C, Coulouris G, Avagyan V, et al. (2009) BLAST+: Architecture and applications. BMC Bioinform 10: 421. https://doi.org/10.1186/1471-2105-10-421
    [37] Cold Spring Harbor Protocols: Lysis buffer for protein extraction (2024). Available from: https://cshprotocols.cshlp.org/content/2014/9/pdb.rec081273.full
    [38] Kulak NA, Pichler G, Paron I, et al. (2014) Minimal, encapsulated proteomic-sample processing applied to copy-number estimation in eukaryotic cells. Nat Methods 11: 319-324. https://doi.org/10.1038/nmeth.2834
    [39] Kovalchuk SI, Jensen ON, Rogowska-Wrzesinska A (2019) FlashPack: fast and simple preparation of ultrahigh-performance capillary columns for LC-MS*[S]. Mol Cell Proteom 18: 383-390. https://doi.org/10.1074/mcp.TIR118.000953
    [40] Tyanova S, Temu T, Cox J (2016) The MaxQuant computational platform for mass spectrometry-based shotgun proteomics. Nat Protoc 11: 2301-2319. https://doi.org/10.1038/nprot.2016.136
    [41] Tyanova S, Temu T, Sinitcyn P, et al. (2016) The Perseus computational platform for comprehensive analysis of (prote) omics data. Nat Methods 13: 731-740. https://doi.org/10.1038/nmeth.3901
    [42] The UniProt Consortium.UniProt: the Universal Protein Knowledgebase in 2023. Nucl Acids Res (2022) 51: D523-D531. https://doi.org/10.1093/nar/gkac1052
    [43] STRING: Protein-Protein Interaction Networks, Functional Enrichment Analysis (2024). Available from: https://string-db.org/cgi/input?sessionId=bJepuA5XwxF1&input_page_show_search=on
    [44] NCBI Basic Local Alignment Search Tool: Protein BLAST (2024). Available from: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins
    [45] Benjamini Y, Hochberg (1995) Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc B: Stat Methodol 57: 289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
    [46] Benjamini Y, Krieger AM, Yekutieli D (2006) Adaptive linear step-up procedures that control the false discovery rate. Biometrika 93: 491-507. https://doi.org/10.1093/biomet/93.3.491
    [47] Staudenmaier H, Van Hove B, Yaraghi Z, et al. (1989) Nucleotide sequences of the fecBCDE genes and locations of the proteins suggest a periplasmic-binding-protein-dependent transport mechanism for iron(III) dicitrate in Escherichia coli. J Bacteriol 171: 2626-2633. https://doi.org/10.1128/jb.171.5.2626-2633.1989
    [48] Brugger SD, Bomar L, Lemon KP (2016) Commensal–pathogen interactions along the human nasal passages. PLoS Pathog 12: e1005633. https://doi.org/10.1371/journal.ppat.1005633
    [49] Penberthy WT, Sadri M, Zempleni J (2020) Biotin. Present knowledge in nutrition. Cambridge: Academic Press 289-303. https://doi.org/10.1016/B978-0-323-66162-1.00017-2
    [50] Hu T, Wei Z, Ju Q, et al. (2021) Sex hormones and acne: State of the art. JDDG 19: 509-515. https://doi.org/10.1111/ddg.14426
    [51] Ray RR, Lahiri D, Chatterjee A, et al. (2021) Bacteria and biofilms as natural inhabitants of our body. Biofilm-Mediated Diseases: Causes and Controls. New York: Springer 47-71. https://doi.org/10.1007/978-981-16-0745-5_3
    [52] Ponomarenko EA, Krasnov GS, Kiseleva OI, et al. (2023) Workability of mRNA sequencing for predicting protein abundance. Genes 14: 2065. https://doi.org/10.3390/genes14112065
    [53] Bakker R, Ellers J, Roelofs D, et al. (2023) Combining time-resolved transcriptomics and proteomics data for Adverse Outcome Pathway refinement in ecotoxicology. Sci Total Environ 869: 161740. https://doi.org/10.1016/j.scitotenv.2023.161740
    [54] Sharma VK, Akavaram S, Bayles DO (2022) Genomewide transcriptional response of Escherichia coli O157:H7 to norepinephrine. BMC Genomics 23: 107. https://doi.org/10.1186/s12864-021-08167-z
    [55] Achermann Y, Tran B, Kang M, et al. (2015) Immunoproteomic identification of in vivo-produced Propionibacterium acnes proteins in a rabbit biofilm infection model. Clin Vaccine Immunol 22: 467-476. https://doi.org/10.1128/CVI.00760-14
    [56] Perraud Q, Kuhn L, Fritsch S, et al. (2020) Opportunistic use of catecholamine neurotransmitters as siderophores to access iron by Pseudomonas aeruginosa. Environ Microbiol 24: 878-893. https://doi.org/10.1111/1462-2920.15372
    [57] Perez-Riverol Y, Bai J, Bandla C, et al. (2022) The PRIDE database resources in 2022: A Hub for mass spectrometry-based proteomics evidences. Nucleic Acids Res 50: D543-D552. https://doi.org/10.1093/nar/gkab1038
    [58] Deutsch EW, Bandeira N, Perez-Riverol Y, et al. (2023) The ProteomeXchange Consortium at 10 years: 2023 update. Nucleic Acids Res 51: D1539-D1548. https://doi.org/10.1093/nar/gkac1040
  • microbiol-10-02-019-s001.pdf
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

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

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

Metrics

Article views(1501) PDF downloads(88) Cited by(0)

Figures and Tables

Figures(4)  /  Tables(10)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog