Pharmacokinetic modeling of Gadolinium nanoparticles (Gd-NPs) with the sojourn time in vasa vasorum for the contrast enhanced MRI

Deposition of lipid in the artery wall called atherosclerosis is recognized as a major cause of cardiovascular disease that leads to death worldwide. A better understanding into factors that may influence the delivery of gadolinium nanoparticles (Gd-NPs) that enhances quality of magnetic resonance imaging in diagnosis may provide a vital key for atherosclerotic treatment. In this study, we propose a delay differential model for describing the dynamics of Gd-NPs in bloodstream, peripheral arteries, and vasa vasorum with two phenomena of Gd-NPs during a sojourn in vasa vasorum. We then investigate the dynamical behaviors of Gd-NPs and explore the effects of sojourn time and transfer rates of Gd-NPs on the concentration of Gd-NPs in vasa vasorum at the 12th hour after the administration of gadolinium chelates contrast media and also the maximum concentration of Gd-NPs in peripheral arteries and vasa vasorum. Our results suggest that the sojourn of Gd-NPs in vasa vasorum may lead to complex behaviors of Gd-NPs dynamics, and transfer rates of Gd-NPs may have a significant impact on the concentration of Gd-NPs.


Introduction
Atherosclerosis, which is a common cardiovascular disease involving inflammation and accumulation of lipid plaques within the arterial wall, is recognized as a major cause of death worldwide [1,2]. There is recent evidence suggesting that vasa vasorum, which is a network of microvessels contained in tunicae adventitia that supply oxygen and nutrients to the arterial wall, is associated with the occurrence of atherosclerotic plaque progression [3,4]. For normal people, the oxygen is usually delivered by vasa vasorum and lumen to respectively reach outer and inner layers of the arterial wall. However, for atherosclerotic patients, vasa vasorum has changed in shape due to angiogenic expansion resulting in it becoming an only supply source for the entire wall [5,6]. Consequently, observation of vasa vasorum expansion can clinically help predict the pathogenesis of atherosclerosis [7].
Magnetic resonance imaging (MRI) is an important technique currently used to diagnose and detect the pathogenesis of early plaque. It is known that contrast agents are clinically employed to help enhance the quality of medical imaging modalities [8,9]. Usually, the contrast agents are presented in forms of nanoparticles (NPs) or liposomes containing metal ions such as iron(II), iron(III), gadolinium(II), gadolinium(III), manganese(II), or manganese(III) [10][11][12]. Gadolinium(III) complexes (gadolinium-based contrast agents (GBCAs) or gadolinium nanoparticles (Gd-NPs)) are the most widely clinically used agent through the intravenous route [13]. The working principle of Gd-NPs to assist the image acquisition is binding between the Gd-NPs and the specific biological environment such as molecular markers of disease. For atherosclerosis, there are α ν β 3 -integrins acting as the receptors of inflammatory endothelial cells on vasa vasorum, and Gd-NPs are designed and synthesized in a variety of chemical methods with the aim of targeting the α ν β 3 -integrin receptors [6,14,15].
Even though the integrin binding helps enhance MR images, there are concerns about the safety of the contrast agents. In order to investigate factors that may influence their delivery and to predict possible effects inside the organs, physiologically-based pharmacokinetic (PBPK) models have often been chosen as important tools in many research studies. To investigate the multi-modality imaging technologies, Bartlett et al. [16] developed a three-compartmental model for plasma, interstitial and intracellular fluids to investigate the biodistribution of transferrin (tf )-targeted and nontargeted siRNA nanoparticles and the cellular uptake of a tumor. The study indicated that the tf-targeted NPs were more effective in attaching to the cell target and can help reduce tumor activity by 50 percent. Barboriak et al. [17] constructed three PBPK models including well-mixed, delay, and dispersion models to study the concentration-time curves (CTCs) of the gadolinium contrast agent at three position organs, which are cerebral cortex, superior sagittal sinus, and psoas muscle. Their results suggest that the best model to predict Gd-CTCs is the dispersion model. Neubauer et al. [18] examined the behavior of perfluorocarbon NPs while they were transported in blood, peripheral, and vasa vasorum by using a three-compartmental pharmacokinetic model. In their work, the NPs have the ability to target α ν β 3 -integrin receptors of atherosclerotic tissue resulting in the enhancement of contrast MR images. The purpose of their study was to develop models that are compatible with the new drug carriers. Wenger et al. [19] proposed a mathematical model based on a mass balance equation to evaluate distribution, accumulation, and excretion of the synthesized two NPs, polyacrylamide (PAA) and polyethylene glycol-coated polyacrylamide (PAApeg). Taheri et al. [20] studied the pharmacokinetics of Gd-DTPA administered at a bolus dose to optimize an area under the curve (AUC) of Gd-DTPA and found that reduction of Gd-DTPA together with optimized acquisition of DCE-MRI results in shorter acquisition time and less exposure of subjects.
In this present study, we incorporate the time-delay terms that reflect the sojourn of Gd-NPs in vasa vasorum into the mathematical model developed [18]. We also explore two types of functions describing circulation of Gd-NPs during the temporary stay in vasa vasorum: with and without natural clearance of Gd-NPs. Our goals are to (i) investigate the dynamical behaviors of Gd-NPs including periodic, asymptomatic, and oscillating solutions in bloodstream, peripheral arteries, and vasa vasorum; and (ii) explore the effects of certain parameters such as the sojourn time that Gd-NPs stay in vasa vasorum and transfer rates between the compartments on the concentration of Gd-NPs in vasa vasorum at the 12th hour after administration and the maximum concentration of Gd-NPs within peripheral arteries and vasa vasorum.

Pharmacokinetic model
Previously, a three-compartmental pharmacokinetic model was proposed to investigate the dynamical behaviors of Gd-NPs involved in the enhancement of MRI of early atherosclerosis [18]. Those three compartments represent bloodstream (y 1 ), peripheral arteries (y 2 ), and vasa vasorum (y 3 ). In this work, the previous model is further extended to incorporate a fixed length of the sojourn time (τ ) that Gd-NPs remain trapped in the third compartment so that the MR image signal remains elevated in comparison to the blood signal. Transfer of Gd-NPs between the first two compartments occurs at rates k 12 and k 21 . Transfer of Gd-NPs between the third compartment and the first compartment occurs at rates k 13 and k 31 . In addition to the distribution phase of Gd-NPs, the model also includes the clearance phase since Gd-NPs are normally removed from the bloodstream by phagocytosis in the reticuloendothelial system. An elimination constant for the clearance phase is denoted by k e . As specific binding of α v β 3 -targeted Gd-NPs to α v β 3 -integrin can take place in vasa vasorum to enhance the MR signal, such transfer may include both passive and active transports of Gd-NPs. Consequently, the dynamics of Gd-NPs in those three compartments can be described by where y 1 (t), y 2 (t), and y 3 (t) represent the concentration of Gd-NPs at time t (μM) in bloodstream, peripheral arteries, and vasa vasorum, respectively. In the model, f is a function of τ where τ represents the sojourn time in the vasa vasorum compartment. Two types of f are considered here to describe two different behaviors of Gd-NPs. The first case is f (τ ) = 1, which reflects the completely retained amount of Gd-NPs during their temporary stay in the vasa vasorum compartment. The second case is f (τ ) = e -k 31 τ , which reflects the decreasing amount of Gd-NPs at rate k 31 during the sojourn stay in the vasa vasorum compartment. To tackle the delay system of equations in (1) numerically, the corresponding non-delayed system is initially solved with the initial conditions where D is the gadolinium dosage of 4.6 μmol Gd/kg body weight injected to each test subject and v 1 is the distribution volume of bloodstream. Such a model is solved until the time reaches t = τ to obtain a new set of initial conditions for (1) on [0, τ ] as follows: where Y 0 i (θ ) for i = 1, 2, 3 denotes a solution of the corresponding non-delayed system at time θ . For t ∈ (τ , ∞), the delay system is solved with conditions in (2). Parameters used in Table 1 Model parameters [18] Name Description Targeted Unit k 12 Transfer rate from peripheral arteries to bloodstream 1.220 1/h k 21 Transfer rate from bloodstream to peripheral arteries 1.290 1/h k 13 Transfer rate from bloodstream to vasa vasorum 1.100 1/h k 31 Transfer rate from vasa vasorum to bloodstream 1.400 1/h k e Clearance rate from bloodstream 0.126 1/h v 1 Volume of blood 0.061 L/kg  Table 1, while the flow diagram for describing the dynamical behaviors of Gd-NPs is illustrated in Fig. 1.
Note that model (1) can be extended to incorporate diffusion and chemotaxis of Gd-NPs. Examples of preceding studies that take into account those terms are such as Li et al. [21] and Viglialoro and Woolley [22].

Analytic results
Some analytic results of the model including an equilibrium and certain stability conditions are demonstrated in this section. Some other works related to the analysis of oscillatory and periodic solutions in detail can be found in [23] and [24]. By setting each equation in (1) equaling zero, the Gd-NPs-free equilibrium is obtained, E 0 = (y 0 1 , y 0 2 , y 0 3 ) = (0, 0, 0).

Theorem 1 The Gd-NPs-free equilibrium is locally asymptotically stable.
Proof First, we compute the eigenvalues (λ) of the Jacobian matrix from the equation, det (J -λI) = 0. The characteristic equation of J is given by with a 1 = k 21 + k 12 + k 13 + k e + k 31 , a 2 = k 21 k 13 + k 21 k e + k 21 k 31 + k 31 k 12 + k 31 k e , It is clear that a 1 > 0, a 3 > 0, and a 1 a 2 > a 3 . Hence, according to the Routh-Hurwitz criteria [25], all eigenvalues of J have negative real parts. Consequently, E 0 is locally asymptotically stable.

τ = 0 and f (τ ) = 1
For this case, as an equilibrium solution does not depend on time, there exists the Gd-NPs-free equilibrium of (1), E 0 . To investigate the stability of E 0 , let us define Then the linearized system of (1) in a matrix form can be described by where A 1 and A 2 are given by The characteristic equation of the linearized system is given by Since τ is transcendental and the system has infinitely many eigenvalues, investigating the stability of E 0 is much more difficult. Hence, we shall instead explore the distribution of roots of Eq. (3).
By Rouche's theorem and the continuity in τ , Eq. (3) has positive real parts if and only if it has purely imaginary roots. Then we study whether we can find purely imaginary roots.
By separating real and imaginary parts and squaring their equations, we obtain Next, let Consequently, Eq. (5) can be rewritten as Proposition 1 Because γ < 0, Eq. (6) has at least one positive root, denoted by Z 0 .

τ = 0 and f (τ ) = e -k 31 τ
Since there is still only a Gd-NPs-free equilibrium (E 0 ) for this case, we follow a similar approach in the previous part to explore stability of the equilibrium. Hence, the charac-teristic equation of the linearized system is given by with b n for n = 1, . . . , 5.
Substituting λ(τ ) = iω(τ ) with ω > 0 into Eq. (7), we have By separating real and imaginary parts in Eq. (8), the real part satisfies the following equation: while the imaginary part satisfies the equation Next, Eq. (10) can be rewritten as Then, substituting Eq. (11) into Eq. (9) yields Notice that Eq. (12)  · e -k 31 τ . Based on parameter ranges of this study, there is no intersection between those two curves. Hence, we cannot find ω that satisfies Eq. (12) (see Fig. 2(a)-(d)). Consequently, there is no periodic solution of the system. Theorem 3 If Eq. (12) has no ω, then E 0 is always locally asymptotically stable for all τ > 0.

Numerical results
In this section, numerical results for dynamical behaviors of Gd-NPs in bloodstream, peripheral arteries, and vasa vasorum are illustrated to support theoretical results obtained from the previous section and to demonstrate the effects of transfer rates and the sojourn time that Gd-NPs remain trapped in the vasa vasorum compartment. Note that the numerical simulations are based on the Runge-Kutta method of order 4 incorporated with the initial conditions for solving delay differential equations [26].
In case f (τ ) = 1 for τ > 0 or equivalently Gd-NPs are trapped in the vasa vasorum at a fixed length of time before escaping to other parts, dynamical behaviors of Gd-NPs according to different values of τ (0.3, 1.7, 1.96, 2.5) are shown in Fig. 3. As a critical value of τ  Table 1 is approximately 1.96 hours, Gd-NPs dynamics change around such a value. If τ < 1.96, the system solution tends to the Gd-NPs-free equilibrium, and it starts oscillating when τ becomes closer to τ 0 (see Fig. 3(a)-(b)). The results correspond to Theorem 2 that the Gd-NPs-free equilibrium of system (1) is locally asymptotically stable when τ < 1.96. At τ = 1.96, Fig. 3(c) shows that the solution fluctuates in a periodic pattern. Those dynamical behaviors are in agreement with Theorem 2 so that a Hopf bifurcation occurs at τ = τ 0 = 1.96. In Fig. 3(d), when τ > 1.96, the solution fluctuates with an increasing amplitude and becomes divergent. This is in agreement with the unstable behaviors around the equilibrium after τ > τ 0 in Theorem 2. In Fig. 4, we demonstrate phase portraits at different values of τ in correspondence with results in   Fig. 4(c). When τ = 2.5, the solution trajectory diverges to infinity over time (see Fig. 4(d)).
Next, we investigate the effects of transfer rates on the concentration of Gd-NPs at the 12th hour after the initial injection according to different values of τ (see Fig. 5). Results in Fig. 5(a) suggest that increasing the transfer rate of Gd-NPs from bloodstream to peripheral arteries may result in the decreased concentration of Gd-NPs in vasa vasorum. However, such a trend of results may slightly be affected by the increasing value of τ . In Fig. 5(b), the transfer rate from peripheral arteries to bloodstream has a significant impact on the concentration of Gd-NPs when it is quite small. When it becomes larger, the concentration gradually increases. The increasing value of τ may alter the trend by delaying the increasing concentration of Gd-NPs. Results in Fig. 5(c) demonstrate the increasing concentration of Gd-NPs according to the transfer rate from bloodstream to vasa vasorum. Figure 5(d) suggests that increasing the transfer rate from vasa vasorum to blood- Moreover, we explore the effects of transfer rates on the maximum concentration of Gd-NPs inside peripheral arteries and vasa vasorum when Gd-NPs do not and do remain in vasa vasorum at the fixed length of time τ = 0, 1.7. Note that based on the preceding study from the preceding work, Gd-NPs remain trapped in vasa vasorum for approximately τ = 1.7 hours. In Fig. 6, the overall results suggest similar trends of results according to transfer rates for both τ = 0 and τ = 1.7, but the concentration of Gd-NPs for the τ = 0 case is generally higher. Our results in Fig. 6(a) demonstrate that when the transfer rate from bloodstream to peripheral arteries increases, the concentration of Gd-NPs in peripheral arteries considerably increases, while the concentration of Gd-NPs in vasa vasorum moderately decreases. Figure 6(b) shows the significantly decreasing concentration of Gd-NPs in peripheral arteries and the slightly increasing concentration of Gd-NPs in vasa vasorum according to the increasing transfer rate from peripheral arteries to bloodstream. In Fig. 6(c), when the transfer rate from bloodstream to vasa vasorum increases, it may result in the increase of Gd-NPs concentration in vasa vasorum but in the moderate decrease of Gd-NPs concentration in peripheral arteries. The increasing transfer rate from vasa vasorum to bloodstream may lead to the decreasing concentration of Gd-NPs in vasa vasorum but the slightly increasing concentration of Gd-NPs in peripheral arteries for τ = 0 and the slightly increasing then decreasing concentration of Gd-NPs in peripheral arteries for τ = 1.7 (see Fig. 6(d)). In addition, when the sojourn time is present, it may also result in the increasing concentration after the reduction of Gd-NPs in vasa vasorum. Our results in Fig. 6(e) demonstrate that the concentration of Gd-NPs in vasa vasorum and peripheral arteries decreases if the elimination rate increases.
In the case f (τ ) = e -k 31 τ for τ = 0 or equivalently Gd-NPs are trapped in the vasa vasorum at a fixed length of time and eliminated during the stay at a rate of e -k 31 τ before escaping to other parts, time-series plots and phase portraits that represent the dynamical behaviors near the Gd-NPs-free equilibrium for τ = 0.3 and τ = 2.5 are illustrated in   Figure 7(a)-(b) demonstrate the concentration of Gd-NPs according to time that Gd-NPs are eventually eliminated and a phase portrait that shows a solution approaching the origin for τ = 0.3. In Fig. 7(c)-(d), our results show that the concentration of Gd-NPs approaches zero according to time and a corresponding phase portrait that shows a solution approaching the origin for τ = 2.5.
We further investigate the effects of transfer rates on the concentration of Gd-NPs at the 12th hour after the initial injection when f (τ ) = e -k 31 τ for τ = 0. Overall, our results suggest that increasing τ leads to the higher concentration of Gd-NPs. In Fig. 8(a), the transfer rate of Gd-NPs from bloodstream to peripheral arteries rarely affects the concentration of Gd-NPs. Figure 8(b) shows that when the transfer rate from peripheral arteries to bloodstream increases, the concentration of Gd-NPs increases at the beginning before becoming saturated according to the transfer rate. The increasing transfer rate of Gd-NPs Figure 9 The maximum concentration of Gd-NPs in peripheral arteries and vasa vasorum without and with particles remaining trapped in vasa vasorum (τ = 0, 1.7) for f (τ ) = e -k 31 τ from bloodstream to vasa vasorum significantly influences the concentration of Gd-NPs and results in higher concentration (see Fig. 8(c)). In contrast, when the transfer rate from vasa vasorum to bloodstream increases, the concentration of Gd-NPs significantly decreases as shown in Fig. 8(d). Figure 8(e) illustrates that the increasing clearance rate may reduce the concentration of Gd-NPs in vasa vasorum. Figure 9(a)-(e) illustrates that the longer sojourn time that Gd-NPs remain trapped generally results in the higher concentration of Gd-NPs in vasa vasorum and the lower concentration of Gd-NPs in peripheral arteries. In Fig. 9(a), increasing the transfer rate of Gd-NPs from bloodstream to peripheral arteries may moderately increase the concentration of Gd-NPs in peripheral arteries and slightly reduce the concentration of Gd-NPs in vasa vasorum. Results in Fig. 9(b) suggest that increasing the transfer rate of Gd-NPs from peripheral arteries to bloodstream may result in the significant decrease of Gd-NPs concentration in peripheral arteries and the small increase of Gd-NPs concentration in vasa vasorum. As illustrated in Fig. 9(c), when the transfer rate of Gd-NPs from bloodstream to vasa vasorum increases, the concentration of Gd-NPs in peripheral arteries intermediately decreases, while the concentration of Gd-NPs in vasa vasorum considerably increases. In Fig. 9(d), increasing the transfer rate of Gd-NPs from vasa vasorum to bloodstream may result in the slightly decreasing concentration of Gd-NPs in peripheral arteries and the significantly increasing concentration of Gd-NPs in vasa vasorum when the sojourn of Gd-NPs in vasa vasorum is present. The increasing elimination rate of Gd-NPs generally leads to the decreasing concentration of Gd-NPs in both peripheral arteries and vasa vasorum (see Fig. 9(e)). When results in Figs. 6 and 9 are compared, the difference of f (τ ) may affect the concentration trends in vasa vasorum. When Gd-NPs remain trapped in vasa vasorum, the concentration of Gd-NPs is lower than when they do not for f (τ ) = 1, while it is higher for f (τ ) = e -k 31 τ .

Conclusion
In this study, the three-compartmental pharmacokinetic model for describing the distribution of Gd-NPs used in the contrast enhanced MRI is extended to incorporate a sojourn period of particle circulation in vasa vasorum. The particle dynamics are then described by a system of delay differential equations. Two types of circulation are considered: with and without natural clearance during the particle sojourn.
When the sojourn of Gd-NPs in vasa vasorum is not included (τ = 0), we show that the particle-free equilibrium point is locally asymptotically stable. When the fixed period of circulation is present (τ = 0), the delay system is analyzed to determine the conditions for asymptomatic stability and the critical time of a periodic solution. Our findings suggest that only when Gd-NPs are not eliminated during the sojourn in vasa vasorum, a periodic solution occurs under certain conditions. If Gd-NPs are eliminated by natural clearance during their stay in vasa vasorum, a system solution always tends to the particle-free equilibrium point.
In addition, we demonstrate the dynamics of Gd-NPs numerically according to the stability conditions. Our results are in agreement with the conditions and show asymptotic, periodic, and unstable solution behaviors of Gd-NPs when there is no clearance of Gd-NPs during their stay in vasa vasorum (f (τ ) = 1). Generally, a solution starts to oscillate when τ approaches its critical value (τ = 1.96), shows a periodic behavior when τ = 1.96, and then shows an unstable behavior when τ > 1.96. Consequently, to ensure that Gd-NPs are completely eliminated, Gd-NPs should be designed not to linger in vasa vasorum for too long. However, if Gd-NPs are naturally cleared during their stay in vasa vasorum (f (τ ) = e -k 31 τ ), the particle-free equilibrium point is asymptotically stable. In such a case, a solution will tend to the equilibrium point or Gd-NPs will always be eventually eliminated from the body. As τ is normally less than the critical value in reality, we further investigate how transfer rates of Gd-NPs affect their presence in blood stream, peripheral arteries, and vasa vasorum after twelve hours of injection, and also the maximum concentration of Gd-NPs. It is found that increasing transfer rates (k 12 , k 31 , and k e ) lowers the concentration of Gd-NPs, while increasing k 21 and k 13 leads to the higher concentration of Gd-NPs. For the maximum concentration, our results suggest that transfer rates affect the concentration of Gd-NPs in peripheral arteries and vasa vasorum. Note that results are similar for both types of Gd-NPs circulation. In addition, overall results suggest that increasing τ may result in the higher concentration of Gd-NPs. Moreover, the maximum of Gd-NPs in vasa vasorum for f (τ ) = 1 is generally lower than that for f (τ ) = e -k 31 τ .
Finally, we believe that this study may help gain a better understanding into the effects of the fixed period of particle circulation and transfer parameters, may also help design a therapy that enhances the efficacy of MRI, and may suggest possible administration strategies that could reduce toxicity from overdose.