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

Toxicity associated with gadolinium-based contrast-enhanced examinations

  • Received: 20 February 2021 Accepted: 26 April 2021 Published: 08 May 2021
  • This article reports known and emerging adverse health effects associated with the administration of gadolinium-based contrast agents. It focuses on the issue of the incomplete excretion of these drugs leading to the deposition of gadolinium in the tissues of the patients. The evidence of deposition is reviewed. The analysis presents gaps in our knowledge but also suggests neglected or still poorly considered parameters to possibly explain discrepancies among studies (e.g. off-label use; rate of administration; gadolinium concentration in the pharmaceutical formulation, cumulative metal toxicity). The article also presents a critical assessment of some aspects reported in the literature as well as future needs. Potential biases in the investigation and evaluation of the health/clinical implications associated with gadolinium deposition are pointed out. The analysis emphasizes that the vast majority of the clinical studies conducted up to date on gadolinium-based contrast agents were designed to assess acute toxicity and diagnostic efficacy of the agents, not to identify long-term health effects.

    Citation: Silvia Maria Lattanzio. Toxicity associated with gadolinium-based contrast-enhanced examinations[J]. AIMS Biophysics, 2021, 8(2): 198-220. doi: 10.3934/biophy.2021015

    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
  • This article reports known and emerging adverse health effects associated with the administration of gadolinium-based contrast agents. It focuses on the issue of the incomplete excretion of these drugs leading to the deposition of gadolinium in the tissues of the patients. The evidence of deposition is reviewed. The analysis presents gaps in our knowledge but also suggests neglected or still poorly considered parameters to possibly explain discrepancies among studies (e.g. off-label use; rate of administration; gadolinium concentration in the pharmaceutical formulation, cumulative metal toxicity). The article also presents a critical assessment of some aspects reported in the literature as well as future needs. Potential biases in the investigation and evaluation of the health/clinical implications associated with gadolinium deposition are pointed out. The analysis emphasizes that the vast majority of the clinical studies conducted up to date on gadolinium-based contrast agents were designed to assess acute toxicity and diagnostic efficacy of the agents, not to identify long-term health effects.



    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.



    Conflict of interest



    The authors declare no conflict of interest.

    [1] Tamburrini O, Aprile I, Falcone C, et al. (2011) Off-label use of intravascular iodinated organic and MR contrast media. Radiol Med 116: 1-14. doi: 10.1007/s11547-010-0601-5
    [2] Meloni MM, Barton S, Xu L, et al. (2017) Contrast agents for cardiovascular magnetic resonance imaging: an overview. J Mater Chem B 5: 5714-5725. doi: 10.1039/C7TB01241A
    [3] Lauffer RB (1987) Paramagnetic metal complexes as water proton relaxation agents for NMR imaging: theory and design. Chem Rev 87: 901-927. doi: 10.1021/cr00081a003
    [4] Schörner W, Kazner E, Laniado M, et al. (1984) Magnetic resonance tomography (MRT) of intracranial tumours: Initial experience with the use of the contrast medium gadolinium-DTPA. Neurosurg Rev 7: 303-312. doi: 10.1007/BF01892910
    [5] Essig M, Anzalone N, Combs SE, et al. (2012) MR imaging of neoplastic central nervous system lesions: review and recommendations for current practice. Am J Neuroradiol 33: 803-817. doi: 10.3174/ajnr.A2640
    [6] Anzalone N, Gerevini S, Scotti R, et al. (2009) Detection of cerebral metastases on magnetic resonance imaging: intraindividual comparison of gadobutrol with gadopentetate dimeglumine. Acta Radiol 50: 933-940. doi: 10.1080/02841850903095385
    [7] Malayeri AA, Brooks KM, Bryant LH, et al. (2016) National Institutes of health perspective on reports of gadolinium deposition in the brain. J Am Coll Radiol 13: 237-241. doi: 10.1016/j.jacr.2015.11.009
    [8] Weinmann HJ, Brasch RC, Press WR, et al. (1984) Characteristics of gadolinium-DTPA complex: a potential NMR contrast agent. Am J Roentgenol 142: 619-624. doi: 10.2214/ajr.142.3.619
    [9] Lin SP, Brown JJ (2007) MR contrast agents: Physical and pharmacologic basics. J Magn Reson Imaging 25: 884-899. doi: 10.1002/jmri.20955
    [10] Ersoy H, Rybicki FJ (2007) Biochemical safety profiles of gadolinium-based extracellular contrast agents and nephrogenic systemic fibrosis. J Magn Reson Imaging 26: 1190-1197. doi: 10.1002/jmri.21135
    [11] Bellin MF, Van Der Molen AJ (2008) Extracellular gadolinium-based contrast media: an overview. Eur J Radiol 66: 160-167. doi: 10.1016/j.ejrad.2008.01.023
    [12] Evans CH (1990)  Biochemistry of the Lanthanides New York: Plenum Press. doi: 10.1007/978-1-4684-8748-0
    [13] De León-Rodríguez LM, Martins AF, Pinho MC, et al. (2015) Basic MR relaxation mechanisms and contrast agent design. J Magn Reson Imaging 42: 545-565. doi: 10.1002/jmri.24787
    [14] Levine D, McDonald RJ, Kressel HY (2018) Gadolinium retention after contrast-enhanced MRI. JAMA 320: 1853-1854. doi: 10.1001/jama.2018.13362
    [15] Shepherd M, Lata S, Mani S, et al. (2009) Anaphylaxis to gadolinium radiocontrast: a case report and review of the literature. J La State Med Soc 161: 282-284.
    [16] Raisch DW, Garg V, Arabyat R, et al. (2014) Anaphylaxis associated with gadolinium-based contrast agents: Data from the food and drug administration's adverse event reporting system and review of case reports in the literature. Expert Opin Drug Saf 13: 15-23. doi: 10.1517/14740338.2013.832752
    [17] Franckenberg S, Berger F, Schaerli S, et al. (2018) Fatal anaphylactic reaction to intravenous gadobutrol, a gadolinium-based MRI contrast agent. Radiol Case Rep 13: 299-301. doi: 10.1016/j.radcr.2017.09.012
    [18] Jung JW, Kangand HR, Kim MH, et al. (2012) Immediate hypersensitivity reaction to gadolinium-based MR contrast media. Radiology 264: 414-422. doi: 10.1148/radiol.12112025
    [19] Morzycki A, Bhatia A, Murphy KJ (2017) Adverse reactions to contrast material: a Canadian update. Can Assoc Radiol J 68: 187-193. doi: 10.1016/j.carj.2016.05.006
    [20] Behzadi AH, Farooq Z, Newhouse JH, et al. (2018) MRI and CT contrast media extravasation a systematic review. Medicine 97: e0055. doi: 10.1097/MD.0000000000010055
    [21] Varela DC, Sepulveda P, Prieto J, et al. (2015) Extravasation of intravenous contrast media: What every radiologist should know. Rev Chil Radiol 21: 151-157. doi: 10.4067/S0717-93082015000400006
    [22] Blasco-Perrin H, Glaser B, Pienkowski M, et al. (2013) Gadolinium induced recurrent acute pancreatitis. Pancreatology 13: 88-89. doi: 10.1016/j.pan.2012.12.002
    [23] Unal O, Arslan H (1999) Cardiac arrest caused by IV gadopentetate dimeglumine. Am J Roentgenol 172: 1141. doi: 10.2214/ajr.172.4.10587169
    [24] US Food and Drug AdministrationMedical Imaging Drugs Advisory Committee Meeting, FDA briefing document: Gadolinium retention after gadolinium based contrast magnetic resonance imaging in patients with normal renal function. (2017) .
    [25] Maramattom BV, Manno EM, Wijdicks EFM, et al. (2005) Gadolinium encephalopathy in a patient with renal failure. Neurology 64: 1276-1278. doi: 10.1212/01.WNL.0000156805.45547.6E
    [26] Hui FK, Mullins M (2009) Persistence of gadolinium contrast enhancement in CSF: A possible harbinger of gadolinium neurotoxicity? Am J Neuroradiol 30: e1. doi: 10.3174/ajnr.A1205
    [27] Kim SH, Jo EJ, Kim MY, et al. (2013) Clinical value of radiocontrast media skin tests as a prescreening and diagnostic tool in hypersensitivity reactions. Ann Allergy Asthma Immunol 110: 258-262. doi: 10.1016/j.anai.2013.01.004
    [28] Shellock FG, Kanal E (1999) Safety of magnetic resonance imaging contrast agents. J Magn Res Im 10: 477-484. doi: 10.1002/(SICI)1522-2586(199909)10:3<477::AID-JMRI33>3.0.CO;2-E
    [29] Grobner T (2006) Gadolinium-a specific trigger for the development of nephrogenic fibrosing dermopathy and nephrogenic systemic fibrosis? Nephrol Dial Transpl 21: 1104-1108. doi: 10.1093/ndt/gfk062
    [30] Thomson LK, Thomson PC, Kingsmore DB, et al. (2015) Diagnosing nephrogenic systemic fibrosis in the post-FDA restriction era. J Magn Reson Imaging 41: 1268-1271. doi: 10.1002/jmri.24664
    [31] Larson KN, Gagnon AL, Darling MD, et al. (2015) Nephrogenic systemic fibrosis manifesting a decade after exposure to gadolinium. JAMA Dermatol 151: 1117-1120. doi: 10.1001/jamadermatol.2015.0976
    [32] Bernstein EJ, Schmidt-Lauber C, Kay J (2012) Nephrogenic systemic fibrosis: a systemic fibrosing disease resulting from gadolinium exposure. Best Pract Res Cl Rheumatol 26: 489-503. doi: 10.1016/j.berh.2012.07.008
    [33] Sanyal S, Marckmann P, Scherer S, et al. (2011) Multiorgan gadolinium (Gd) deposition and fibrosis ina patient with nephrogenic systemic fibrosis-an autopsy-based review. Nephrol Dial Transplant 26: 3616-3626. doi: 10.1093/ndt/gfr085
    [34] Kay J, Bazari H, Avery LL, et al. (2008) Case 6-2008: a 46-year-old woman with renal failure and stiffness of the joints and skin. New Engl J Med 358: 827-838. doi: 10.1056/NEJMcpc0708697
    [35] Bhave G, Lewis JB, Chang SS (2008) Association of gadolinium based magnetic resonance imaging contrast agents and nephrogenic systemic fibrosis. J Urol 180: 830-835. doi: 10.1016/j.juro.2008.05.005
    [36] Cowper SE, Su LD, Bhawan J, et al. (2001) Nephrogenic fibrosing dermopathy. Am J Dermatopath 23: 383-393. doi: 10.1097/00000372-200110000-00001
    [37] Cowper SE, Robin HS, Steinberg SM, et al. (2000) Scleromyxedema-like cutaneous disease in renal dialysis patients. Lancet 356: 1000-1001. doi: 10.1016/S0140-6736(00)02694-5
    [38] Weigle JP, Broome DR (2008) Nephrogenic systemic fibrosis: chronic imaging findings and review of the medical literature. Skeletal Radiol 37: 457-464. doi: 10.1007/s00256-008-0464-1
    [39] Zou Z, Zhang HL, Roditi GH, et al. (2011) Nephrogenic systemic fibrosis: review of 370 biopsy-confirmed cases. JACC: Cardiovasc Imag 4: 1206-1216. doi: 10.1016/j.jcmg.2011.08.013
    [40] Morris MF, Zhang Y, Zhang H, et al. (2009) Features of nephrogenic systemic fibrosis on radiology examinations. Am J Roentgenol 193: 61-69. doi: 10.2214/AJR.08.1352
    [41] Tsushima Y, Kanal E, Thomsen HS (2010) Nephrogenic systemic fibrosis: risk factors suggested from Japanese published cases. Brit J Radiol 83: 590-595. doi: 10.1259/bjr/17689538
    [42] Thomsen HS, Morcos SK, Almén T, et al. (2013) Nephrogenic systemic fibrosis and gadolinium-based contrast media: updated ESUR contrast medium safety committee guidelines. Eur Radiol 23: 307-318. doi: 10.1007/s00330-012-2597-9
    [43] Mazhar SM, Shiehmorteza M, Kohl CA, et al. (2009) Nephrogenic systemic fibrosis in liver disease: a systematic review. J Magn Reson Imaging 30: 1313-1322. doi: 10.1002/jmri.21983
    [44] Elmholdt TR, Jørgensen B, Ramsing M, et al. (2010) Two cases of nephrogenic systemic fibrosis after exposure to the macrocyclic compound gadobutrol. NDT Plus 3: 285-287.
    [45] Kay J (2008) Gadolinium and nephrogenic systemic fibrosis: The evidence of things not seen. Clev Clin J Med 75: 112. doi: 10.3949/ccjm.75.2.112
    [46] Todd DJ, Kay J (2008) Nephrogenic systemic fibrosis: an epidemic of gadolinium toxicity. Curr Rheumatol Rep 10: 195-204. doi: 10.1007/s11926-008-0033-6
    [47] Sieber MA, Lengsfeld P, Frenzel T, et al. (2008) Preclinical investigation to compare different gadolinium-based contrast agents regarding their propensity to release gadolinium in vivo and to trigger nephrogenic systemic fibrosis-like lesions. Eur Radiol 18: 2164-2173. doi: 10.1007/s00330-008-0977-y
    [48] Semelka RC, Prybylskib JP, Ramalho M (2019) Influence of excess ligand on nephrogenic systemic fibrosis associated with nonionic, linear gadolinium-based contrast agents. Magn Res Imaging 58: 174-178. doi: 10.1016/j.mri.2018.11.015
    [49] US Food and Drug Administration FDA request boxes warning for contrast agents used to improve MRI images (2007) .Available from: http://wayback.archive-it.org/7993/20170112033008/http://www.fda.gov/NewsEvents/Newsroom/PressAnnouncements/2007/ucm108919.htm.
    [50] Khawaja AZ, Cassidy DB, Al Shakarchi J, et al. (2015) Revisiting the risks of MRI with Gadolinium based contrast agents: review of literature and guidelines. Insights Imaging 6: 553-558. doi: 10.1007/s13244-015-0420-2
    [51] Canga A, Kislikova M, Martínez-Gálvez M, et al. (2014) Renal function, nephrogenic systemic fibrosis and other adverse reactions associated with gadolinium-based contrast media. Nefrologia 34: 428-438.
    [52] Martin DR, Krishnamoorthy SK, Kalb B, et al. (2010) Decreased incidence of NSF in patients on dialysis after changing gadolinium contrast-enhanced MRI protocols. J Magn Reson Imaging 31: 440-446. doi: 10.1002/jmri.22024
    [53] Altun E, Martin DR, Wertman R, et al. (2009) Nephrogenic systemic fibrosis: Change in incidence following a switch in gadolinium agents and adoption of a gadolinium policy—report from two US universities. Radiology 253: 689-696. doi: 10.1148/radiol.2533090649
    [54] Xia D, Davis RL, Crawford JA, et al. (2010) Gadolinium released from MR contrast agents is deposited in brain tumors: in situ demonstration using scanning electron microscopy with energy dispersive X-ray spectroscopy. Acta Radiol 51: 1126-1136. doi: 10.3109/02841851.2010.515614
    [55] Kanda T, Ishii K, Kawaguchi H, et al. (2014) High signal intensity in the dentate nucleus and globus pallidus on unenhanced T1-weighted MR images: Relationship with increasing cumulative dose of a gadolinium based contrast material. Radiology 270: 834-841. doi: 10.1148/radiol.13131669
    [56] McDonald RJ, McDonald JS, Kallmes DF, et al. (2015) Intracranial gadolinium deposition after contrast-enhanced MR imaging. Radiology 275: 772-782. doi: 10.1148/radiol.15150025
    [57] Olchowy C, Cebulski K, Łasecki M, et al. (2017) The presence of the gadolinium-based contrast agent depositions in the brain and symptoms of gadolinium neurotoxicity-a systematic review. PLoS One 12: e0171704. doi: 10.1371/journal.pone.0171704
    [58] Pullicino R, Radon M, Biswas S, et al. (2018) A review of the current evidence on gadolinium deposition in the brain. Clin Neuroradiol 28: 159-169. doi: 10.1007/s00062-018-0678-0
    [59] Gianolio E, Gregorio ED, Aime S (2019) Chemical insights into the issues of Gd retention in the brain andother tissues upon the administration of Gd-containing MRI contrast agents. Eur J Inorg Chem 2019: 137-151. doi: 10.1002/ejic.201801220
    [60] Kanda T, Osawa M, Oba H, et al. (2015) High signal intensity in dentate nucleus on unenhanced T1-weightedMR images: association with linear versus macrocyclic gadolinium chelate administration. Radiology 275: 803-809. doi: 10.1148/radiol.14140364
    [61] Weberling LD, Kieslich PJ, Kickingereder P, et al. (2015) Increased signal intensity in the dentate nucleus on unenhanced T1-weighted images after gadobenate dimeglumine administration. Invest Radiol 50: 743-748. doi: 10.1097/RLI.0000000000000206
    [62] Errante Y, Cirimele V, Mallio CA, et al. (2014) Progressive increase of T1 signal intensity of the dentate nucleus on unenhanced magnetic resonance images is associated with cumulative doses of intravenously administered gadodiamide in patients with normal renal function, suggesting dechelation. Invest Radiol 49: 685-690. doi: 10.1097/RLI.0000000000000072
    [63] Zhang Y, Cao Y, Shih GL, et al. (2017) Extent of signal hyperintensity on unenhanced T1-weighted brain MR images after more than 35 administrations of linear gadolinium-based contrast agents. Radiology 282: 516-525. doi: 10.1148/radiol.2016152864
    [64] Radbruch A, Weberling LD, Kieslich PJ, et al. (2015) Gadolinium retention in the dentate nucleus and globus pallidus is dependent on the class of contrast agent. Radiology 275: 783-791. doi: 10.1148/radiol.2015150337
    [65] Malhotra A, LeSar B, Wu X, et al. (2018) Progressive T1 shortening of the dentate nucleus in patients with multiple sclerosis: Result of multiple administrations of linear gadolinium contrast agents versus intrinsic disease. Am J Roentgenol 211: 1099-1105. doi: 10.2214/AJR.17.19155
    [66] Flood TF, Stence NV, Maloney JA, et al. (2017) Pediatric brain: Repeated exposure to linear gadolinium-based contrast material is associated with increased signal intensity at unenhanced T1-weighted MR imaging. Radiology 282: 222-228. doi: 10.1148/radiol.2016160356
    [67] Adin ME, Kleinberg L, Vaidya D, et al. (2015) Hyperintense dentate nuclei on T1-weighted MRI: relation to repeat gadolinium administration. Am J Neuroradiol 36: 1859-1865. doi: 10.3174/ajnr.A4378
    [68] Miller JH, Hu HH, Pokorney A, et al. (2015) MRI brain signal intensity changes of a child during the course of 35 gadolinium contrast examinations. Pediatrics 136: e1637-e1640. doi: 10.1542/peds.2015-2222
    [69] Hu HH, Pokorney A, Towbin RB, et al. (2016) Increased signal intensities in the dentate nucleus and globus pallidus on unenhanced T1-weighted images: evidence in children undergoing multiple gadolinium MRI exams. Pediatr Radiol 46: 1590-1598. doi: 10.1007/s00247-016-3646-3
    [70] Bae S, Lee H, Han K, et al. (2017) Gadolinium deposition in the brain: association with various GBCAs using a generalized additive model. Eur Radiol 27: 3353-3361. doi: 10.1007/s00330-016-4724-5
    [71] Quattrocchi CC, Mallio CA, Errante Y, et al. (2015) Gadodiamide and dentate nucleus T1 hyperintensity in patients with meningioma evaluated by multiple follow-up contrast-enhanced magnetic resonance examinations with no systemic interval therapy. Invest Radiol 50: 470-472. doi: 10.1097/RLI.0000000000000154
    [72] McDonald JS, McDonald RJ, Jentoft ME, et al. (2017) Intracranial gadolinium deposition following gadodiamide-enhanced magnetic resonance imaging in pediatric patients: a case-control study. JAMA Pediatr 171: 705-707. doi: 10.1001/jamapediatrics.2017.0264
    [73] Mallio CA, Vullo GL, Messina L, et al. (2020) Increased T1 signal intensity of the anterior pituitary gland on unenhanced magnetic resonance images after chronic exposure to gadodiamide. Invest Radiol 55: 25-29. doi: 10.1097/RLI.0000000000000604
    [74] Gianolio E, Bardini P, Arena F, et al. (2017) Gadolinium retention in the rat brain: Assessment of the amounts of insoluble gadolinium-containing species and intact gadolinium complexes after repeated administration of gadolinium-based contrast agents. Radiology 285: 839-849. doi: 10.1148/radiol.2017162857
    [75] Rasschaert M, Schroeder JA, Wu TD, et al. (2018) Multimodal imaging study of gadolinium presence in rat cerebellum: differences between Gd chelates, presence in the Virchow-Robin space, association with lipofuscin, and hypotheses about distribution pathway. Invest Radiol 53: 518. doi: 10.1097/RLI.0000000000000490
    [76] Radbruch A, Richter H, Fingerhu S, et al. (2019) Gadolinium deposition in the brain in a large animal model. Comparison of linear and macrocyclic gadolinium-based contrast agents. Invest Radiol 54: 531-536. doi: 10.1097/RLI.0000000000000575
    [77] Boyken J, Frenzel T, Lohrke J, et al. (2018) Gadolinium accumulation in the deep cerebellar nuclei and globus pallidus after exposure to linear but not macrocyclic gadolinium-based contrast agents in a retrospective pig study with high similarity to clinical conditions. Invest Radiol 53: 278-285. doi: 10.1097/RLI.0000000000000440
    [78] Robert P, Violas X, Grand S, et al. (2016) Linear gadolinium-based contrast agents are associated with brain gadolinium retention in healthy rats. Invest Radiol 51: 73-82. doi: 10.1097/RLI.0000000000000241
    [79] Strzeminska I, Factor C, Robert P, et al. (2020) Long-term evaluation of gadolinium retention in rat brain after single injection of a clinically relevant dose of gadolinium-based contrast agents. Invest Radiol 55: 138-143. doi: 10.1097/RLI.0000000000000623
    [80] Radbruch A, Haase R, Kieslich PJ, et al. (2017) No signal intensity increase in the dentate nucleus on unenhanced T1-weighted MR images after more than 20 serial injections of macrocyclic gadolinium-based contrast agents. Radiology 282: 699-707. doi: 10.1148/radiol.2016162241
    [81] Radbruch A, Haase R, Kickingereder P, et al. (2017) Pediatric brain: no increased signal intensity in the dentate nucleus on unenhanced T1-weighted MR images after consecutive exposure to a macrocyclic gadolinium-based contrast agent. Radiology 283: 828-836. doi: 10.1148/radiol.2017162980
    [82] Tibussek D, Rademacher C, Caspers J, et al. (2007) Gadolinium brain deposition after macrocyclic gadolinium administration: a pediatric case-control study. Radiology 285: 223-230. doi: 10.1148/radiol.2017161151
    [83] Schneider GK, Stroeder J, Roditi G, et al. (2017) T1 signal measurements in pediatric brain: findings after multiple exposures to gadobenate dimeglumine for imaging of non neurologic disease. Am J Neuroradiol 38: 1799-1806. doi: 10.3174/ajnr.A5270
    [84] Conte G, Preda L, Cocorocchio E, et al. (2017) Signal intensity change on unenhanced T1-weighted images in dentate nucleus and globus pallidus after multiple administrations of gadoxetate disodium: an intraindividual comparative study. Eur Radiol 27: 4372-4378. doi: 10.1007/s00330-017-4810-3
    [85] Ryu YJ, Choi YH, Cheon J, et al. (2018) Pediatric brain: Gadolinium deposition in dentate nucleus and globus pallidus on unenhanced T1-weighted images is dependent on the type of contrast agent. Invest Radiol 53: 246-255. doi: 10.1097/RLI.0000000000000436
    [86] Stanescu AL, Shaw DW, Murata N, et al. (2020) Brain tissue gadolinium retention in pediatric patients after contrast-enhanced magnetic resonance exams: pathological confirmation. Pediatr Radiol 50: 388-396. doi: 10.1007/s00247-019-04535-w
    [87] Bjørnerud A, Vatnehol SAS, Larsson C, et al. (2017) Signal enhancement of the dentate nucleus at unenhanced MR imaging after very high cumulative doses of the macrocyclic gadolinium-based contrast agent gadobutrol: an observational study. Radiology 285: 434-444. doi: 10.1148/radiol.2017170391
    [88] Splendiani A, Perri M, Marsecano C, et al. (2018) Effects of serial macrocyclic based contrast materials gadoterate meglumine and gadobutrol administrations on gadolinium related dentate nuclei signal increases in unenhanced t1-weighted brain: a retrospective study in 158 multiple sclerosis (MS) patients. Radiol Med 123: 125-134. doi: 10.1007/s11547-017-0816-9
    [89] Stojanov DA, Aracki-Trenkic A, Vojinovic S, et al. (2016) Increasing signal intensity within the dentate nucleus and globus pallidus on unenhanced T1W magnetic resonance images in patients with relapsing remitting multiple sclerosis: Correlation with cumulative dose of a macrocyclic gadolinium-based contrast agent, gadobutrol. Eur Radiol 26: 807-815. doi: 10.1007/s00330-015-3879-9
    [90] Tedeschi E, Palma G, Canna A, et al. (2016) In vivo dentate nucleus MRI relaxometry correlates with previous administration of gadolinium-based contrast agents. Eur Radiol 26: 4577-4584. doi: 10.1007/s00330-016-4245-2
    [91] Lattanzio SM, Imbesi F (2020) Fibromyalgia associated with repeated gadolinium contrast-enhanced MRI examinations. Radiol Case Rep 15: 534-541. doi: 10.1016/j.radcr.2020.02.002
    [92] Roberts DR, Welsh CA, LeBel DP, et al. (2017) Distribution map of gadolinium deposition within the cerebellum following GBCA administration. Neurology 88: 1206-1208. doi: 10.1212/WNL.0000000000003735
    [93] Gibby WA, Gibby KA, Gibby WA (2004) Comparison of Gd-DTPA-BMA (Omniscan) versus Gd-HPDO3A (ProHance) retention in human bone tissue by inductively coupled plasma atomic emission spectroscopy. Invest Radiol 39: 138-142. doi: 10.1097/01.rli.0000112789.57341.01
    [94] White GW, Gibby WA, Tweedle MF (2006) Comparison of Gd(DTPA-BMA)(Omniscan) versus Gd(HPDO3A)(ProHance) relative to gadolinium retention in human bone tissue by inductively coupled mass spectroscopy. Invest Radiol 41: 272-278. doi: 10.1097/01.rli.0000186569.32408.95
    [95] Darrah TH, Prutsman-Pfeiffer JJ, Poreda RJ, et al. (2009) Incorporation of excess gadolinium into human bone from medical contrast agents. Metallomics 1: 479-488. doi: 10.1039/b905145g
    [96] Murata N, Gonzalez-Cuyar LF, Murata K, et al. (2016) Macrocyclic and other non–group 1 gadolinium contrast agents deposit low levels of gadolinium in brain and bone tissue: Preliminary results from 9 patients with normal renal function. Invest Radiol 51: 447-53. doi: 10.1097/RLI.0000000000000252
    [97] Lord ML, Chettle DR, Gräfe JL, et al. (2018) Observed deposition of gadolinium in bone using a new noninvasive in vivo biomedical device: Results of a small pilot feasibility study. Radiology 287: 96-103. doi: 10.1148/radiol.2017171161
    [98] Turyanskaya A, Rauwol M, Pichler V, et al. (2020) Detection and imaging of gadolinium accumulation inhuman bone tissue by micro- and submicro-XRF. Sci Rep 10: 6301. doi: 10.1038/s41598-020-63325-9
    [99] Vidaud C, Bourgeois D, Meyer D (2012) Bone as target organ for metals: the case of f-elements. Chem Res Toxicol 25: 1161-1175. doi: 10.1021/tx300064m
    [100] Gräfe JL, McNeill FE (2018) Measurement of gadolinium retention: current status and review from an applied radiation physics perspective. Physiol Meas 39: 06TR01. doi: 10.1088/1361-6579/aacc16
    [101] Hasegawa M, Duncan BR, Marshall DA, et al. (2020) Human hair as a possible surrogate marker of retained tissue gadolinium. A pilot autopsy study correlating gadolinium concentrations in hair with brain and other tissues among decedents who received gadolinium-based contrast agents. Invest Radiol 55: 636-642. doi: 10.1097/RLI.0000000000000681
    [102] Saussereau E, Lacroix C, Cattaneo A, et al. (2008) Hair and fingernail gadolinium ICP-MS contents in an overdose case associated with nephrogenic systemic fibrosis. Forensic Sci Int 176: 54-57. doi: 10.1016/j.forsciint.2007.06.026
    [103]  US Food and Drug Administration, 5-18, 2015 Available from: https://www.fda.gov/drugs/drug-safety-and-availability/fda-drug-safety-communication-fda-evaluating-risk-brain-deposits-repeated-use-gadolinium-based.
    [104]  EMA/625317/2017. EMA's final opinion confirms restrictions on use of linear gadolinium agents in body scans Available at: http://www.ema.europa.eu/docs/en _ GB/document _ library/ Referrals _ document/gadolinium _ contrast _ agents _ 31/ European _ Commission _ final _ decision/WC500240575.pdf.
    [105] Lancelot E, Desché P (2020) Gadolinium retention as a safety signal: experience of a manufacturer. Invest Radiol 55: 20-24. doi: 10.1097/RLI.0000000000000605
    [106]  PMDA, Revision of Precautions, Gadodiamide hydrate Meglumine gadopentetate, 2017 Available from: http://www.pmda.go.jp/files/000221377.pdf.
    [107]  PMDA, Revision of Precautions, Gadoxetate sodium, Gadoteridol, Meglumine gadoterate, Gadobutrol, 2017 Available from: http://www.pmda.go.jp/files/000221376.pdf.
    [108] Kanda T, Nakai Y, Hagiwara A, et al. (2017) Distribution and chemical forms of gadolinium in the brain: a review. Br J Radiol 90: 20170115. doi: 10.1259/bjr.20170115
    [109] Bracco Diagnostics Bayer, Guerbet GE Healthcare Important drug warning for all gadolinium-based contrast agents [dear health care provider letter] (2018) .Available from: https://www.guerbet.com/media/uh4h4kon/dhcp-letter-05-02-2018-signed.pdf.
    [110] Harvey HB, Gowda V, Cheng G (2019) Gadolinium deposition disease: a new risk management threat. J Am Coll Radiol 17: 546-550. doi: 10.1016/j.jacr.2019.11.009
    [111] Semelka RC, Commander CW, Jay M, et al. (2016) Presumed gadolinium toxicity in subjects with normal renal function a report of 4 cases. Invest Radiol 51: 661-665. doi: 10.1097/RLI.0000000000000318
    [112] Semelka RC, Ramalho M, Jay M (2016) Summary of special issue on gadolinium bioeffects and toxicity with a look to the future. Magn Reson Imaging 34: 1399-1401. doi: 10.1016/j.mri.2016.09.002
    [113] US Food and Drug AdministrationMedical Imaging Drugs Advisory Committee Meeting, FDA briefing document: Gadolinium retention after gadolinium-based contrast magnetic resonance imaging in patients with normal renal function, 27–28. (2017) .
    [114] Burke LMB, Ramalho M, AlObaidy M, et al. (2016) Self-reported gadolinium toxicity: A survey of patients with chronic symptoms. Magn Reson Imaging 34: 1078-1080. doi: 10.1016/j.mri.2016.05.005
    [115] Semelka RC, Ramalho J, Vakharia A, et al. (2016) Gadolinium deposition disease: initial description of adisease that has been around for a while: a family of disorders. Magn Res Imaging 34: 1383-1390. doi: 10.1016/j.mri.2016.07.016
    [116] Roberts DR, Lindhorst SM, Welsh CT, et al. (2016) High levels of gadolinium deposition in the skin of a patient with normal renal function. Invest Radiol 51: 280-289. doi: 10.1097/RLI.0000000000000266
    [117] Barbieri S, Schroeder C, Froehlich JM, et al. (2016) High signal intensity in dentate nucleus and globus pallidus on unenhanced T1-weighted MR images in three patients with impaired renal function and vascular calcification. Contrast Media Mol Imaging 11: 245-250. doi: 10.1002/cmmi.1683
    [118] Swaminathan S (2016) Gadolinium toxicity: iron and ferroportin as central targets. Magnet Reson Imaging 34: 1373-1376. doi: 10.1016/j.mri.2016.08.016
    [119] Di Gregorio ED, Furlan C, Atlante S, et al. (2020) Gadolinium retention in erithrocytes and leukocytes from human and murine blood upon treatment with gadolinium-based contrast agents for magnetic resonance imaging. Invest Radiol 55: 30-37. doi: 10.1097/RLI.0000000000000608
    [120] Di Gregorio E, Ferrauto G, Furlan C, et al. (2018) The issue of gadolinium retained in tissue. Invest Radiol 53: 167-172. doi: 10.1097/RLI.0000000000000423
    [121] Kartamihardja AAP, Hanaoka H, Andriana P, et al. (2019) Quantitative analysis of Gd in the protein content of the brain following single injection of gadolinium-based contrast agents (GBCAs) by size exclusion chromatography. Br J Radiol 92: 20190062. doi: 10.1259/bjr.20190062
    [122] Newton BB, Jimenez SA (2009) Mechanism of NSF: New evidence challenging the prevailing theory. J Magn Reson Imaging 30: 1277-1283. doi: 10.1002/jmri.21980
    [123] Taoka T, Jost G, Frenzel T, et al. (2018) Impact of the glymphatic system on the kinetic and distribution of gadodiamide in rat brain. Invest Radiol 53: 529-534. doi: 10.1097/RLI.0000000000000473
    [124] Nehra AK, McDonald RJ, Bluhm AM, et al. (2018) Accumulation of gadolinium in human cerebrospinal fluid after gadobutrol-enhanced MR imaging: a prospective observational cohort study. Radiology 288: 416-423. doi: 10.1148/radiol.2018171105
    [125] McDonald RJ, Levine D, Weinreb J, et al. (2018) Gadolinium retention: A research roadmap from the 2018 NIH/ACR/RSNA workshop on gadolinium chelates. Radiology 289: 517-534. doi: 10.1148/radiol.2018181151
    [126] Le Fur M, Caravan P (2019) The biological fate of gadolinium-based MRI contrast agents: a call to action for bioinorganic chemists. Metallomics 11: 240-254. doi: 10.1039/C8MT00302E
    [127] Tweedle MF (2016) Gadolinium deposition: Is it chelated or dissociated gadolinium? How can we tell? Magn Res Imaging 34: 1377-1382. doi: 10.1016/j.mri.2016.09.003
    [128] Kiviniemi A, Gardberg M, Ek P, et al. (2019) Gadolinium retention in gliomas and adjacent normal brain tissue: association with tumor contrast enhancement and linear/macrocyclic agents. Neuroradiology 61: 535-544. doi: 10.1007/s00234-019-02172-6
    [129] Kanda T, Fukusato T, Matsuda M, et al. (2015) Gadolinium-based contrast agent accumulates in the brain even in subjects without severe renal dysfunction: evaluation of autopsy brain specimens with inductively coupled plasma mass spectroscopy. Radiology 276: 228-232. doi: 10.1148/radiol.2015142690
    [130] Herculano-Houzel S (2009) The human brain in numbers: a linearly scaled-up primate brain. Front Hum Neurosci 3: 31. doi: 10.3389/neuro.09.031.2009
    [131] Popescu BFG, Robinson CA, Rajput A, et al. (2009) Iron, copper, and zinc distribution of the cerebellum. Cerebellum 8: 74-79. doi: 10.1007/s12311-008-0091-3
    [132] Kromrey ML, Liedtke KR, Ittermann T, et al. (2017) Intravenous injection of gadobutrol in an epidemiological study group did not lead to a difference in relative signal intensities of certain brain structures after 5years. Eur Radiol 27: 772-777. doi: 10.1007/s00330-016-4418-z
    [133] Staks T, Schuhmann-Giampieri G, Frenzel T, et al. (1994) Pharmacokinetics, dose proportionality and tolerability of gadobutrol after single intravenous injection in healthy volunteers. Invest Radiol 29: 709-715. doi: 10.1097/00004424-199407000-00008
    [134] Gutierrez JE, Rosenberg M, Duhaney M, et al. (2015) Phase 3 efficacy and safety trial of gadobutrol, a 1.0 molar macrocyclic MR imaging contrast agent, in patients referred for contrast-enhanced MR imagingof the central nervous system. J Magn Reson Imaging 41: 788-796. doi: 10.1002/jmri.24583
    [135] Kuwatsuru R, Takahashi S, Umeoka S, et al. (2015) A multicenter, randomized, controlled, single-blind comparison phase III study to determine the efficacy and safety of gadobutrol 1.0 M versus gadopentetate dimeglumine following single injection in patients referred for contrast-enhanced MRI of the body regions or extremities. J Magn Reson Imaging 41: 404-413. doi: 10.1002/jmri.24566
    [136] Liang Z, Ma L, Wang D, et al. (2012) Efficacy and safety of gadobutrol (1.0 M) versus gadopentetate dimeglumine (0.5 M) for enhanced MRI of CNS lesions: A phase III, multicenter, single-blind, randomized study in Chinese patients. Mag Res Insights 5: MRI-S9348.
    [137] Naito S, Tazaki H, Okamoto T, et al. (2017) Comparison of nephrotoxicity between two gadolinium-contrasts, gadodiamide and gadopentetate in patients with mildly diminished renal failure. J Toxicol Sci 42: 379-384. doi: 10.2131/jts.42.379
    [138] Semelka RC, Hernandes MA, Stallings CG, et al. (2013) Objective evaluation of acute adverse events andimage quality of gadolinium-based contrast agents (gadobutrol and gadobenate dimeglumine) by blinded evaluation. Pilot study. Magn Reson Imaging 31: 96-101. doi: 10.1016/j.mri.2012.06.025
    [139] Tanaka A, Masumoto T, Yamada H, et al. (2016) A Japanese, multicenter, open-label, phase 3 study to investigate the safety and efficacy of gadobutrol for contrast-enhanced MR imaging of the central nervous system. Magn Reson Med Sci 15: 227-236. doi: 10.2463/mrms.mp.2015-0083
    [140] Zech CJ, Schwenke C, Endrikat J (2019) Diagnostic efficacy and safety of gadoxetate disodium vs gadobenate dimeglumine in patients with known or suspected focal liver lesions: Results of a clinical phase III study. Magn Reson Insight 12: 1178623X19827976.
    [141] Tweedle MF, Wedeking P, Kumar K (1995) Biodistribution of radiolabeled, formulated gadopentetate, gadoteridol, gadoterate, and gadodiamide in mice and rats. Invest Radiol 30: 372-380. doi: 10.1097/00004424-199506000-00008
    [142] Rocklage SM, Worah D, Kim SH (1991) Metal ion release from paramagnetic chelates: what is tolerable? Magn Reson Med 22: 216-221. doi: 10.1002/mrm.1910220211
    [143] Khairinisa MA, Takatsuru Y, Amano I, et al. (2018) The effect of perinatal gadolinium-based contrast agents. Invest Radiol 53: 110-118. doi: 10.1097/RLI.0000000000000417
    [144] Ray JG, Vermeulen MJ, Bharatha A, et al. (2016) Association between MRI exposure during pregnancy and fetal and childhood outcomes. JAMA 316: 952-961. doi: 10.1001/jama.2016.12126
    [145] Runge VM, Kuehl TJ, Jackson CB, et al. (2005) Subchronic toxicity of the gadolinium chelates. Acad Radiol 12: S6-S9. doi: 10.1016/j.acra.2005.02.015
    [146] Alkhunizi SM, Fakhoury M, Abou-Kheir W, et al. (2020) Gadolinium retention in the central and peripheral nervous system: implications for pain, cognition, and neurogenesis. Radiology 297: 407-416. doi: 10.1148/radiol.2020192645
    [147] Wang S, Hesse B, Roman M, et al. (2019) Increased retention of gadolinium in the inflamed brain after repeated administration of gadopentetate dimeglumine. Invest Radiol 54: 617-626. doi: 10.1097/RLI.0000000000000571
    [148] Reimer P, Vosshenrich R (2008) Off-label use of contrast agents. Eur Radiol 18: 1096-1101. doi: 10.1007/s00330-008-0886-0
    [149] Essig M, Shiroishi MS, Nguyen TB, et al. (2013) Perfusion MRI: The five most frequently asked technical questions. Am J Roentgenol 200: 24-34. doi: 10.2214/AJR.12.9543
    [150] Wolansky LJ, Cadavid D, Punia V, et al. (2015) Hypophosphatemia is associated with the serial administration of triple-dose gadolinium to patients for brain MRI. J Neuroimaging 25: 379-383. doi: 10.1111/jon.12241
    [151] Essig M, Giesel E, Le-Huu M, et al. (2004) Perfusion MRI in CNS disease: current concepts. Neuroradiology 46: S201-S207. doi: 10.1007/s00234-004-1331-y
    [152] Lee JY, Park JE, Kim HS, et al. (2017) Up to 52 administrations of macrocyclic ionic MR contrast agent are not associated with intracranial gadolinium deposition: multifactorial analysis in 385 patients. PloS One 12: e0183916. doi: 10.1371/journal.pone.0183916
    [153] Ng KH, Ahmad AC, Nizam M, et al.Magnetic resonance imaging: Health effects and safety, Proceedings of the international conference on non-ionizing radiation at UNITEN, Electromagnetic Fields and our Health. (2003) .
    [154] Cho S, Lee Y, Choi YJ, et al. (2014) Enhanced cytotoxic and genotoxic effects of gadolinium following ELF-EMF irradiation in human lymphocytes. Drug Chem Toxicol 37: 440-447. doi: 10.3109/01480545.2013.879662
    [155] Sadiq S, Ghazala Z, Chowdhury A, et al. (2012) Metal toxicity at the synapse: presynaptic, postsynaptic, and long-term effects. J Toxicol 2012: 132671. doi: 10.1155/2012/132671
    [156]  Food U S, Drug Administration, Safety Announcement 2017 Available from: https://www.fda.gov/media/109825/download.
    [157] Veiga M, Mattiazzi P, de Goisc JS, et al. (2020) Presence of other rare earth metals in gadolinium-based contrast agents. Talanta 216: 120940. doi: 10.1016/j.talanta.2020.120940
    [158] Parfrey P (2005) The clinical epidemiology of contrast-induced nephropathy. Cardiovasc Intervent Radiol 28: S3-S11. doi: 10.1007/s00270-005-0196-8
    [159] Karcaaltincaba M, Oguz B, Haliloglu M (2009) Current status of contrast-induced nephropathy and nephrogenic systemic fibrosis in children. Pediatr Radiol 39: S382-S384. doi: 10.1007/s00247-009-1236-3
    [160] Kulaksiz S, Bau M (2011) Anthropogenic gadolinium as a micro-contaminant in tap water used as drinking water in urban areas and megacities. Appl Geochem 26: 1877-1885. doi: 10.1016/j.apgeochem.2011.06.011
    [161] Hatje V, Bruland KW, Flegal AR (2016) Increases in anthropogenic gadolinium anomalies and rare earth element concentrations in San Francisco bay over a 20 year record. Environ Sci Technol 50: 4159-4168. doi: 10.1021/acs.est.5b04322
    [162] Rabiet M, Brissaud F, Seidel JL, et al. (2009) Positive gadolinium anomalies in wastewater treatment plant effluents and aquatic environment in the Hérault watershed (South France). Chemosphere 75: 1057-1064. doi: 10.1016/j.chemosphere.2009.01.036
    [163] Chen Y, Cao XD, Lu Y, et al. (2000) Effects of rare earth metal ions and their EDTA complexes on antioxidant enzymes of fish liver. B Environ Contam Tox 65: 357-365. doi: 10.1007/s001280000136
    [164] Henriques B, Coppola F, Monteiro R, et al. (2019) Toxicological assessment of anthropogenic gadolinium in seawater: Biochemical effects in mussels mytilus galloprovincialis. Sci Total Environ 664: 626-634. doi: 10.1016/j.scitotenv.2019.01.341
    [165] Hanana H, Turcotte P, André C, et al. (2017) Comparative study of the effects of gadolinium chloride and gadolinium-based magnetic resonance imaging contrast agent on freshwater mussel, dreissena polymorphaChemosphere 181: 197-207. doi: 10.1016/j.chemosphere.2017.04.073
    [166] Martino C, Costa C, Roccheria MC, et al. (2018) Gadolinium perturbs expression of skeletogenic genes, calcium uptake and larval development in phylogenetically distant sea urchin species. Aquat Toxicol 194: 57-66. doi: 10.1016/j.aquatox.2017.11.004
    [167] Schmidt K, Bau M, Merschel G, et al. (2019) Anthropogenic gadolinium in tap water and in tap-based beverages from fast-food franchises in six major cities in Germany. Sci Total Environ 687: 1401-1408. doi: 10.1016/j.scitotenv.2019.07.075
  • 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(9913) PDF downloads(483) Cited by(7)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog