Numerical Simulation of Thrombolysis in Robot-Assisted Retinal Vein Cannulation

Beijing Key Laboratory of Advanced Manufacturing Technology, Beijing University of Technology, Beijing 100124, China Faculty of Information Engineering and Automation, Kunming University of Science and Technology, Kunming 650031, China College of Future Technology, University of Chinese Academy of Sciences, Beijing 100049, China Department of Pathology, Sichuan Cancer Hospital & Institute, Sichuan Cancer Center, School of Medicine, University of Electronic Science and Technology of China, Chengdu 610041, China


Introduction
Retinal vein occlusion (RVO) is one of the most common retinal disease where thrombus inside retinal veins distort the delivery of oxygen, resulting in central retinal vein occlusion (CRVO) and branch retinal vein occlusion (BRVO) [1][2][3][4][5]. It impairs the vision of an estimated 16.4 million people in the world [1]. Up to present, the therapeutic strategy for retinal vein occlusion has aimed to limit the damage caused by the obstruction, but cannot dissolve the occlusion in vessel. Consequently, RVO patients need to be treated on a periodic basis with costly therapies to offer limited symptom relief [6,7], such as intravitreal injections [8] and laser photocoagulation [9]. A promising potentially curative treatment method is retinal vein cannulation, also known as retinal endovascular surgery (REVS), as shown in Figure 1. By injecting an anticoagulant directly into the occluded vessel, REVS is focused on thrombolysis to cure RVO completely [10]. Cannulation of the involved vessel and removal of the clot may provide a more permanent solution with a less demanding follow-up [11].
Because of the small size of retinal vessel, surgical robot and machine vision detection [12] are essential for cannulation. Smet et al. injected ocriplasmin into porcine eyes to clear clot by robot-assisted cannulation [13], and Willekens et al. performed cannulation to release drugs for more than 3 min and dissolved clot in 15 of 18 porcine eyes with the help of surgical robot [14]. eir research is limited to animal experiments. However, the shape and size of porcine eyes are different from those of human eyes. Gijbels et al. successfully cured four patients with robot devices as the first robotassisted retinal rein cannulation in the world, but they did not focus on the details during surgery process, such as injection velocity, drug concentration, injection position, and variable size of clot [15].
us, cannulation is still immature and dangerous for smaller vessel. Due to the lack of human operation and experimental data, the environment and important parameters of cannulation remain unclear.
Some models were established to investigate the process of fibrinolysis. Scott L. Diamond et al. modeled clot lysis by solving a set of reaction-diffusion-convection equations which contains no adjustable parameters [16]. Piebalgs et al. developed a multiphysics model for the dissolution of fibrin clots to obtain the shape and subsequent lysis with three different pressure drops [17]. Wootton et al. incorporated the effect of outer convection by solving the continuity and Navier-Stokes equations for blood flow [19]. Recently, stochastic models have been suggested as a method to analyse the effects of clot structure on thrombolysis [20,21]. In addition, Rijken et al. used a purified experiment system to compare the fibrinolytic properties of two molecular forms of tPA [22], and Anand et al. used a numerical model to predict lysis fronts moving across biogels [27]. However, these methods did not consider the critical parameters of cannulation.
is study aims to detect the influence of each parameter, such as injection velocity, drug concentration, and injection position, on the thrombolysis process in robotassisted surgery by theoretical simulation and fluid calculation of finite element. In addition, this study considers the influence of microneedle and variable clot size in this model and cannula into the thrombus when inadvertent operation occurs. e results obtained from this work will be used to optimize the operation process with robot devices in the near future. Figure 2, the computational geometry is a vein with occlusion where the inner diameter is chosen to be 100 μm in accordance with the dimensions of the human retinal vein [23] and the distance from blood vessel inlet to outlet is long enough to reduce the influence of artificial boundary conditions on flow and protein concentration. Additionally, the length of clot is chosen to be from 0.1 mm to 0.6 mm in this simulation. Finally, the implement for cannulation is a glass microneedle that has 500 μm long tip being angulated 30°with respect to the shaft and the diameter of 30 μm (Clunbury Scientific, LLC, BloomfieldHills (MI), USA) [7]. Figure 3 shows the reaction kinetics of thrombolysis that plasminogen (PLG) is activated by tissue-type plasminogen activator (tPA) absorbed on fibrin surface (bound phase) to form plasmin (PLS), which can degrade the fibrin mesh that stabilizes blood clots. α f specie α (tPA, PLS, PLG) in free phase, α b is specie α in bound phase, PLS AP is PLS which cannot dissolve the clot due to inhabitation of AP. Clot consists of a numbers of protofibril which can be degraded by PLS. Bound phase is available to absorb protease on fibrin surface. In order to evaluate the temporal and spatial changes in concentration, this method describes transport equations in the free phase and bound phase.

Reaction Kinetics and Transport Equations of Species.
In order to improving the accuracy of the method, we considered the main inhibitory effect of α2-antiplasmin (α2-AP) in free phase which can inhibit the fibrinolysis extremely [25] and be showed as follows:  x Figure 3: e process of thrombolysis.  Complexity where k AP is antiplasmin reaction coefficient listed in Table 1. Because of the weak effect, this method ignored the influence of α 2 -macroglobulin (MG) and PAI-1 [26]. For the free phase, the equation is written as follows: where (1) term is the rate of change of the free phase, (2) term is the effect of convection on concentration in the free phase, (3) term is the effect of diffusion, (4) term is the rate of absorption in free phase, and (5) term is about the rate of generation of species α (tPA, PLS, PLG). S α and C α are the bound and free phase concentrations of species α, respectively, θ α is the concentration of total species α that can be absorbed to the surface of fibrin, and ε is porosity (voidage) of clot. ese refer to spatial variations. D α,ij is the hydrodynamic dispersion, k ads,α and k rev,α , and K is forward adsorption coefficient, reversible adsorption coefficient, and solubilization constant. ese parameters are invariant constants listed in Table 1.
For the bound phase, the equation is written as follows: where (1) term is the rate of change of concentrations of specie α with effect of convection in bound phase, (2) term is the effect of diffusion, (3) term is the rate of absorption in free phase, and (4) term is about the rate of generation of specie α. e concentration of the total amount of binding sites that are sterically accessible to the proteins can be found by calculating the amount of binding sites per fibrin monomer and the total amount of these monomers in a protofibril: where χ is the number of protofibrils that are sterically accessible, q α is the number of binding sites in a fibrin monomer for species α, and U 0 is volume of computational cell. In this study, we assume that all the protofibrils in the fibrin fiber are sterically accessible: where r 0 is the average radius of a protofibril, which is shown in Table 1. e transport of PLG in the bound phase should subtract the consumption forming PLS. On the contrary, the transport of PLS should plus the generation. e equation is written as follows: where R PLG is the reduction of PLG, G PLS is generation of PLS, and k 2 and k M are Michaelis-Menten reaction coefficients listed in Table 1.

Fibrinolysis of Clot.
In Diamond SL and Anand S's work [16], this method describes the clot as fibrous porous medium composed of fibrin strands where the specific permeability and porosity is correlated with the fiber radius. Clot lysis was described as a homogeneous shrinking of the fiber radius over time [27]: where R f is the fiber radius, L is the concentration of fibrin lysed, and U 0 is volume of computational cell. In addition, R f0 is initial radius of fibrin fiber, N AV is Avogadro's constant, p 0 is amount of protofibril per fiber cross-sectional area, L pf is total length of protofibrils, and α 0 is amount of fibrin monomers per protofibril length, and these constants are shown in Table 1. e cutting of the clot lysed is described as follows: where c is the solubilization rate per number of cuts required by PLS to cleave one unit of fibrin and k cat is lysis reaction rate coefficient. ese constants are shown in Table 1.
In this method, porosity ε can present the decent of thrombolysis and is written as follows: where ε 0 is initial value of voidage, ρ β and ρ fibre is the density of fibrin in clot and fiber, respectively. ese constants are shown in Table 1. In subsequent comparisons, ε is an important dissolution index. e clot lysis is completed when the porosity gets to one. In order to obtain temporal changes in viscous resistance in porous zone, permeability of it needs to be taken: where k is the permeability, n is the permeability coefficient, and λ is the permeability constant.

Governing Equations of Fluid Flow.
In the study, we considered the influence of microneedle and described the blood flow through a porous medium as laminar and heat exchange is ignored [17]. e existence of viscous force leads to alternating vortex which affects flow velocity and pressure 4 Complexity change extremely around microneedle. e time-dependent control equation is where ρ is blood density, μ is viscosity coefficient, u is velocity, p is pressure, and k is permeability of porous medium [18]. We assume the blood flow as Newtonian and incompressible, and clot as isentropic and homogeneous.

Finite Element Analysis.
is work established the multiphysics continuum model to demonstrate species transport and reaction which simulates the entire lytic process of the occlusive clot. e species reaction function was described by equations (1)-(10) and fluent movement was described by equations (11)-(13).
Source terms of tPA, PLG, and PLS implemented in ANSYS Fluent via user-defined functions (UDFs). is study set the blood flow as laminar flow and describes species transport by Species Transport Model where fibrinolysis of clot occurs in porous zone. e multiphysics continuum model was divided into four cell zones of fluid inlet, fluid outlet, needle zone, and clot zone, where subroutine of sources terms was imported as shown in Figure 2. For cell zone conditions, equations (9) and (13) were used to describe fluid porosity and viscous resistance of porous zone in cot zone, respectively. is model set equations (10) and (13), where ε is equal to ε 0 as initialization, and equations (1)- (9) as adjust in hooks of UDFs. In order to show the process of fibrinolysis of clot and the results, the profiles of ε and C α were added into UDFs. e mesh quality of this model is obtained in the CFD simulation, the minimum body size of the global grid is 1.8456e − 3 mm, and the computational grid cells are 274958.
For parameters, we introduce the effect of AP and set the concentration of AP to be 1.0 μM in free phase [16]. For other initial conditions, fibrin density (ρ fibre ) is 0.28 mg/ml [16] and fibrin radius (R f0 ) is 250 nm [28]. e concentration of PLG in the free phase and bound phase is set as 21.65 μM and 34.9 μM, respectively [16]. For the proportion of drug, the concentration of PLG in drug is 2 μM same as concentration of PLG in plasma. And the plasma concentration of tPA is low as 70 nM; thus, the variety of tPA concentration in drug can make influence for dissolution efficiency. is model added these parameters as initialization of hooks of UDFs.
For boundary conditions, the pressure of vessel inlet is 10 Pa and the pressure of vessel outlet is 0 Pa [29]. In order to get the influence of injection velocity, it varies from 10 mm/s to 100 mm/s. In order to simplify the model, nonpermeable and nonslip walls are introduced.
Based on the simulation method which uses transient analysis to calculate the flow at different times, this study sets several simulations to investigate the effect of tPA concentration, injection velocity, injection position, and size of clot.
For tPA concentration, this study sets tPA concentration of the drug between 10 nM and 100 nM to test and obtain the influence of tPA concentration in the drug. As a basic set, the velocity of injection is 40 mm/s, the length of clot is 0.1 mm, and the injection site (the distance from clot to outlet of microneedle) is 0.3 mm in front of the thrombus.
For injection velocity, it changes from 10 mm/s to 100 mm/s. Other conditions remain unchanged, length of clot is 0.1 mm, and the injection site is 0.3 mm in front of the thrombus. However, tPA concentration for injection velocity is constant for high efficiency which is illustrated from the experiment for tPA concentration.
For injection position, it can be divided into external and internal clot. In actual surgery, the position of injection is in front of clot; thus, we set the injection site to 0.1 mm, 0.2 mm, 0.3 mm, 0.4 mm, and 0.5 mm in front of clot, respectively. Additionally, due to the high difficulty and immaturity for cannulation surgery, inadvertent operation may occur which leads cannula into the thrombus. e event of the microneedle being injected inside of the clot was taken into consideration, and three measurements of injection depth along X-axis were used for experiment: 0.

Validation of the Method.
For the sake of validation, this study repeats simulations to replicate the computational and experimental conditions from literatures [22,27]. e results are consistent with those simulated in the literature [27] with k M and k 2 obtained from Table 1. Using the results obtained by Blinc et al. [25], we replicated the experimental conditions from literatures [22] and assumed that the kinetic parameters for plasminogen activation were kM � 2.42 μM and k2 � 0.22 s − 1 [22]. e difference between our results and numerical results of the literature [27] and experimental results of the literature [22] is 13% and 7%, respectively, as shown in Figure 4(a). Figure 4 shows the time taken when the clot is dissolved up to 95% and deviation of comparisons between simulation results and experimental results. Due to a low deviation rate in Figure 4(b), the accuracy of this model can be verified.

Effect of tPA Concentration.
e drug tPA concentration was set as 10 nM to 100 nM to the effects of tPA concentration. Figure 5 shows lysis time over variations of tPA concentration, which explores the influence of tPA e changing trend of nodes is narrow with a sharp peak that rapidly falls off and the profile becomes wider over tPA concentration. We analyzed the data to establish the polynomial fitting curve which can be written as follows: a 1 , a 2 , a 3 , and k 1 and other values can be obtained from Table 2, T tPA is fibrinolytic time, I tPA is tPA concentration of injection, and its R-squared is 0.99784 which means fitting well. is curve is presented in Figure 5 showing a decrease in time for lysis when tPA concentration increased from 0 nM to 50 nM. However, the rate of change of lysis time slows down gradually when tPA concentration reaches over 50 nM. When tPA concentration increased from 0 nM to 50 nM, the 1st derivative of this curve increases rapidly. However, then it increases slowly with concentration reaching over 50 nM, which clearly demonstrates that the most efficient value of tPA concentration is around 50 nM (for thrombolysis of 230 s). Figure 6 shows the fibrinolysis contours of voidage under the condition, where tPA concentration and injection velocity are 50 nM and 40 mm/s, respectively. Figure 6 show the transversion of thrombolysis contours for the 0.1 mm clot where the initial voidage is about 0.98 which means the intact clot and ultimate voidage is 1.0 which means the completed lysis. e time taken for clot lysis to complete with tPA concentration of 50 nM recorded at around 236 s, as shown in Figure 6, where lysis was found to start at the center in front of clot and formed a conical space allowing thrombolysis to undertake gradually. e result matches the simulation data from [21] to further prove accuracy. Figure 7 shows lysis time over various injection velocity with tPA concentration of 50 nM obtained from tPA concentration stimulations. e time profile falls off rapidly and gets wider, and lysis time for clot lysis to complete changed from 273 s to 120 s corresponding to 10 mm/s to 60 mm/s. en, the various scale of dissolution time is reduced from 120 s to 107 s with injection velocity of 60 mm/s to 100 mm/s. e data of simulation and the polynomial fitting curve are presented in Figure 6. And the polynomial fitting equation can be written as follows:

Effect of Injection Velocity.
where T vel is fibrinolytic time over injection velocity, vis injection velocity, its R-squared is 0.99914, which means fitting well, and other parameters are shown in Table 3. e 1st derivative of this curve increased rapidly from to 10 mm/ s to 60 mm/s and slowly from 60 mm/s to 100 mm/s. is profile revealed that the most efficient value of injection velocity is around 60 mm/s (for thrombolysis of 120 s).
In addition, Figure 8 shows  Figures 9 and 10 show the thrombolysis time over the thrombolysis ratio for different positions. In Figure 9, the trend of three profiles is similar where time increases slowly when the thrombolysis ratio is under 0.6, and the increase of time becomes exponential when the thrombolysis ratio is higher than 0.6. e ratio cannot get to 100%, which means that thrombolysis is incomplete when the distance between the microneedle inlet and the front of the clot is under 0.2 mm. In order to explain this phenomenon, Figure 11 shows the tPA concentration contours for the thrombolytic process when the microneedle is 0.2 mm in front of the clot. In Figure 10, the increment of time increases rapidly when the thrombolysis ratio is lower than 0.1; then, the profiles get wider and increase sharply again until thrombolysis ratio reaches 1.0 exploring dissolution of clot completely. rombolysis can be completed with the distance of 0.3 mm and longer distance, and the thrombolysis time is the minimum with the distance of 0.4 mm. Table 4 shows the clot lysis time for internal cannulation. With the increasing length of insertion, the percentage of dissolution drops with shorter time. Figure 12 shows the voidage contours of 0.025 mm inside the clot and 0.1 mm, 0.2 mm, and 0.3 mm in the front of clot to present the evolution of clot lysis. Figure 12 shows clot lysis starts from the port of microneedle and the area of dissolution changes from a cylinder to a funnel to the end of thrombolysis and most of the thrombus around the needle cannot be dissolved. Figures 12(b) and 12(c) show the incomplete thrombolysis in the end, while floccules of clot.

Effect of Clot
Size. Based on the results of condition simulations and injection mode simulations for thrombolysis, this study obtained more efficient values of 50 nM for tPA concentration, 60 mm/s for injection velocity, and 0.4 mm for distance in the front of the clot. In Figure 6, the thrombolysis process is explored by lysis contours for different sizes of 0.1 mm, 0.2 mm, and 0.3 mm. Figure 13 shows thrombolytic time gets longer as tPA concentration gets higher with any clot size between 0.1 mm and 0.6 mm under the condition, where injection velocity is 60 mm/s and injection position is 0.4 mm in front of clot. Figure 14 shows the simulation results of thrombolysis time gets longer over different injection velocities under condition of tPA concentration of 50 nM and injection position of 0.4 mm.

Discussion
In this paper, we obtained the results of four groups of simulations showing the influence of tPA concentration, injection velocity, injection position, and various clot sizes on thrombolysis. From these results, we can explore how the thrombolytic completion and time vary with the critical parameters, and these change rules can be used as a reference for surgery.
Proper tPA concentration is important because of the potential retinal toxicity of high concentration [31]. e curve in Figure 5 can be used to obtain determination rules of drug concentration to be applied to robot-assisted cannulation and is similar to the results in Sriram Anand's work, where they simulated tPA concentration of from 0.1 ng/ml (0.0014 nM) to 100 ng/ml(1.4 nM) with fibrin concentration of 5.88 μm [27]. e reasons of efficient tPA concentration of 50 nM are the strong influence of concentration in the low range (10 nM-50 nM) and decreasing derivative of lysis time changing. is prediction is consistent with the results from Mark W and Johnson MD' work [30], which revealed proper tPA concentration performed well in progressive intraoperative clot dissolution, while higher concentration failed to make efficient thrombolysis.
To investigate thoroughly the effect of inlet velocity, we compared the clot lysis process of 20 mm/s, 50 mm/s, and 80 mm/s by velocity streamlines displaying in Figure 8  20 mm/s causes longer clot lysis time. However, the vessel can be filled with drug at the proper velocity of 50 mm/s, which leads to efficient thrombolysis. is illustrates that appropriately setting the velocity at 60 mm/s is more efficient than having a higher or lower velocity [16]. In Tameesh's study [32], a high injection flow rate of 0.05 ml/min causes   For the injection position, it can be seen that thrombolysis efficiency is good when the thrombolytic ratio is about 0.6 with distance of lower than 0.2 mm. With the increase of thrombolytic ratio, the thrombolytic time increases rapidly, and the thrombolytic process cannot be completed Figure 9. We can get the reason from Figure 11, which appears that the drug cannot penetrate into the thrombus on both sides of the blood vessel wall and flows out with the flow of blood after the clot is dissolved incompletely. e minimum thrombolytic time occurred around the distance of 0.4 mm because there are floccules formed in the end of thrombolytic process with distance of 0.3 mm or longer in Figure 12; thus, the undissolved clot can be washed away by blood flow [33]. Under the current conditions, we can get the position of the highest efficiency is between 0.3 mm and 0.5 mm, which conforms with the results in the literature which illustrates that the best location for cannulation was directly proximal to the occlusion [7].
It can be concluded that the thrombolytic effect is not good when the microneedle is inserted into the thrombus. With low injection velocity, the drug cannot reach the thrombus on the left side and near the vessel wall. Even if the injection velocity increases, we cannot dissolve the thrombus near the blood vessel wall. Because the drug flows out directly with the blood stream after the thrombus on the right side of the micro-needle has been dissolved as Figure 12 shown. Avoiding thrombus insertion, the operation is easier to perform and the blood vessel is avoided to be damaged.
For the simulation results of different sizes of clot with various tPA concentration and injection velocity, we can conclude that, as injection velocity and drug concentration increased, thrombolytic efficiency increased extremely with a specific clot size. When injection velocity or drug concentration goes beyond a range, however, the influence is gradually diminished. e result of tPA concentration are very similar to those described in Hauptmann and Glusa's work [24]. In addition, the trend in change of fibrinolytic time for low injection velocity consists with the results of Piebalgs and Xu's work [17]. Because of this, we can obtain that the higher efficient range of tPA concentration and injection velocity are from 20 nM to 70 nM and from 40 mm/ s to 60 mm/s, respectively. Proper values need to be determined in the simulation database according to the specific thrombus size.
By data analysis, we can get change rules of effects of drug concentration, injection velocity, and injection position on clot lysis. In the process of robot-assisted surgery, Visual Recognition System can obtain the size and position of clot and the environment of injection with deep learning. On the basis of rules and size of clot, optimal tPA concentration and injection velocity can be determined to be fed into the cannulation operating system to increase the efficiency and success of surgery. In addition, cannulation will be performed in a reasonable position in terms of the size of clot to avoid the haemorrhage reaction [7]. is method can guide the catheter tip accurately to the surface of the vein and avoid inserting into clot to perform surgery successfully.
In this study, the simulation method has a number of limitations. is method simplifies the physical model where it ignores the changes of thrombus shape and vessel diameter and idealizes the injection position of the microneedle. Because of this, the convection of blood and drugs may be some deviation from the actual situation, which can affect fibrinolytic time to a certain extent. Constant vessel diameter also leads to limitations of simulation database.
In the future, we plan to supplement the effect of vessel diameter and geometry of clot on thrombolysis to make the simulation database more comprehensive. And these results will be used to assist the performance of cannulation.

Conclusions
Because of few cases of robotic intubation surgery, there is a lack of data on the effects of drug concentration, injection velocity, injection location, and clot size on fibrinolytic efficacy. e accuracy of the multiphysics continuum model established by finite analysis element is verified by comparing with literatures. In this paper, the data are obtained by simulation of cannulation operation where the influence of microneedle is considered. For clot length of 0.1 mm to 0.6 mm, the data indicates that the optimal tPA concentration ranged from 20 nM to 70 nM and the optimal injection velocity ranged from 40 mm/s to 60 mm/s, and the thrombolysis cannot be completed when the microneedle is inserted into the thrombus and the position of the highest efficiency is between 0.3 mm and 0.5 mm under conditions in this paper. In the future, we will consider the influence of vessel diameter to supplement the database and apply the data in this paper to robotic-assisted intubation surgery.

Data Availability
All data generated or used during the study appear in the submitted article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.