Method for simulating diabetic retinopathy process

Abstract

The invention discloses a diabetic retinopathy process simulation method. The method comprises the steps that a two-dimensional discrete coordinate system is constructed, and a blood vessel index function describing unit grid attributes of unit grids is defined; updating the blood glucose concentration, the oxygen concentration, the fibronectin concentration and the growth factor concentration of each unit grid every preset time t; predicting the movement direction of each blood vessel tip in the t time period according to the fibronectin concentration and the growth factor concentration; calculating a first predicted extension distance and a second predicted extension distance according to the motion direction in the current t time period and the motion direction in the last t time period; and updating the unit grid attribute of each unit grid in the two-dimensional discrete coordinate system according to the movement direction in the t time period, the blood vessel existence time of the blood vessel tip and the movement extension distance calculated by the blood vessel tip. According to the method, the oxygen concentration and the blood glucose concentration are incorporated into the model for simulating vascular extension, and accurate simulation prediction of diabetic retinopathy is achieved.

Description

Technical Field

[0001] The invention relates to the technical field of biological simulation, and in particular to a method for simulating a diabetic retinopathy process. Background Art

[0002] Diabetic retinopathy (DR) is one of the common complications of diabetic patients. The retina is the most light-sensitive tissue in the eye. Studies have shown that its health depends on stable blood and oxygen supply. Long-term hyperglycemia can damage retinal blood vessels, cause retinal microcirculation disorders and vascular lesions, lead to retinal hypoxia, hemorrhage, exudation and other problems, destroy the stable blood and oxygen supply environment of the retina, and cause retinopathy. It can be seen that the high blood sugar environment of diabetic patients and the oxygen concentration in the high blood sugar environment are important pathological factors leading to diabetic retinopathy. Establishing a mathematical model that includes the effects of oxygen concentration and blood sugar is of great significance for a deep understanding of the mechanism of diabetic retinopathy and the development of new treatments.

[0003] The simulation of diabetic retinopathy has important scientific significance and clinical application value. The realization of the simulation of diabetic retinopathy process helps to deeply understand the disease mechanism, optimize the treatment plan, discover new treatment methods, and provide support for personalized medicine and medical education.

[0004] Currently, research has been devoted to understanding the pathological mechanism of diabetic retinopathy and developing effective treatments; among them, the role of growth factors (such as VEGF) in retinal angiogenesis has been widely studied. However, for diabetic patients, a high blood sugar environment is an important symptom. In the process of diabetic retinopathy research, the retinal angiogenesis process predicted without considering the blood sugar level of diabetic patients is far from the actual situation. Summary of the invention

[0005] In order to solve the problem that existing studies devoted to understanding the pathological mechanism of diabetic retinopathy and developing effective treatment methods do not consider the blood glucose level of diabetic patients, resulting in a great difference between the predicted results and the actual situation, an embodiment of the present invention provides a method for simulating the diabetic retinopathy process.

[0006] The specific technical solutions are as follows:

[0007] The embodiment of the present invention provides a method for simulating a diabetic retinopathy process, comprising the following steps:

[0008] Step S100, constructing a two-dimensional discrete coordinate system, and defining a blood vessel index function n(x, y, t) that describes and records the unit grid attributes of each unit grid in the two-dimensional discrete coordinate system at each preset time t; wherein the unit grid attributes include blood vessels, blood vessel tips, and environment;

[0009] Step S200, initializing the blood glucose concentration and oxygen concentration of each unit grid in the two-dimensional discrete coordinate system according to the patient's examination results, and initializing the unit grid properties, fibronectin concentration and growth factor concentration of each unit grid in the two-dimensional discrete coordinate system according to a preset initial state;

[0010] Step S300, updating the blood glucose concentration, the oxygen concentration, the fibronectin concentration, and the growth factor concentration of each unit grid at every preset time t according to the current blood glucose concentration of the unit grid, the control equation of the blood glucose concentration, the current oxygen concentration, the control equation of the oxygen concentration, the current unit grid properties, the current fibronectin concentration, the control equation of the fibronectin concentration, the current growth factor concentration, and the growth factor diffusion reaction equation;

[0011] Step S400, predicting the movement direction of each blood vessel tip in the two-dimensional discrete coordinate system within the current time period t according to the fibronectin concentration of the unit grid after t time, the growth factor concentration after t time, and the blood vessel tip movement direction control equation;

[0012] Step S500, predicting the first movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the first movement direction according to the control equation of the fibronectin concentration at the blood vessel tip after t time, the first movement direction and the movement speed of the blood vessel tip, and predicting the second movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the second movement direction according to the control equation of the fibronectin concentration at the blood vessel tip after t time, the second movement direction and the movement speed of the blood vessel tip; wherein the first movement direction is the movement direction of the blood vessel tip in the current t time period of the blood vessel tip; and the second movement direction is the movement direction of the blood vessel tip in the last t time period of the blood vessel tip;

[0013] Step S600: According to the discrete coordinates (x 0 ,y 0 ) and the first movement speed to calculate the first predicted extension distance of each blood vessel tip in the first movement direction in the two-dimensional discrete coordinate system, and the discrete coordinates (x 0 ,y 0) and the second movement speed to calculate the second predicted extension distance of each blood vessel tip in the two-dimensional discrete coordinate system in the second movement direction; wherein the first predicted extension distance or the second predicted extension distance is the distance from the blood vessel tip to the discrete coordinate (x) of the unit grid to which the blood vessel tip belongs in the first movement direction or the second movement direction at the first movement speed or the second movement speed. 0 ,y 0 ) starts moving and reaches the end coordinate (x 1 ,y 1 ) passes through all the unit grids;

[0014] Step S700, performing extension simulation on each blood vessel tip in the two-dimensional discrete coordinate system within this time period t according to the first movement direction of the blood vessel tip, the existence time of the blood vessel tip, the second movement direction of the blood vessel tip, the first movement speed, the second movement speed, the first predicted extension distance and the second predicted extension distance, and updating the unit grid attributes and the blood vessel tip existence time of each unit grid in the two-dimensional discrete coordinate system with the extension simulation of each blood vessel tip, and updating the simulation time when all blood vessel tips in the two-dimensional discrete coordinate system complete the extension simulation; wherein the extension simulation includes blood vessel extension, blood vessel branching occurring during blood vessel extension and blood vessel fusion occurring during blood vessel extension.

[0015] Optionally, the step S300 includes the following steps:

[0016] Predict the blood glucose concentration of each unit grid after time t according to the current blood glucose concentration of the unit grid and the control equation of the blood glucose concentration; and save it under the information of the corresponding unit grid;

[0017] Predict the oxygen concentration of each unit grid after time t based on the current oxygen concentration of the retinal unit grid, the blood glucose concentration after time t, and the control equation of the oxygen concentration; and save it under the information of the corresponding unit grid;

[0018] Predicting the fibronectin concentration of each retinal unit grid after time t according to the current unit grid properties of the retinal unit grid, the current fibronectin concentration and the control equation of the fibronectin concentration, and saving it under the information of the corresponding unit grid;

[0019] According to the current oxygen concentration, growth factor concentration and growth factor diffusion reaction equation of the unit grid, the growth factor concentration of each unit grid after time t is predicted, and the information of the corresponding unit grid is saved.

[0020] Optionally, the control equation of the oxygen concentration is as follows:

[0021]

[0022] in:

[0023] c o is the current oxygen concentration of the retinal cell grid;

[0024] D o is the diffusion coefficient of oxygen in retinal blood vessels, which is a fixed constant;

[0025] F(c o ) is the oxygen consumption rate equation of the retinal unit grid:

[0026]

[0027] A is the normal oxygen consumption rate of the unit grid, which is a fixed constant;

[0028] B is the constant of oxygen consumption of the unit grid under the influence of blood glucose;

[0029] c olim The lowest oxygen concentration for normal physiological activities of the unit grid is a fixed constant. Unit grids below this constant value are in a hypoxic state.

[0030] G(ρ) is the amount of oxygen transported from the surrounding blood vessels in the unit grid. Vascular transport is the source of oxygen for the retinal unit grid. According to the vascular density at the current time, the calculation equation of G(ρ) is as follows:

[0031]

[0032] α o is the oxygen transport rate constant in blood vessels;

[0033] C is the transport efficiency constant under the influence of blood glucose;

[0034] ρ(x,y,t 0 ) is the unit grid (x, y) of the retina at the current time t 0 The blood vessel density at the time is specifically the number of unit grids with blood vessels in the current unit grid (x, y) and the nine unit grids adjacent to the current unit grid.

[0035] Optionally, the control equation of the fibronectin concentration is:

[0036]

[0037] in:

[0038] c f is the current fibronectin concentration of the retinal unit grid; the initial value has been set in step S200;

[0039] ω is the production rate of fibronectin in endothelial cells, which is a fixed constant;

[0040] is the uptake rate of fibronectin by endothelial cells and is a fixed constant.

[0041] Optionally, the growth factor diffusion reaction equation is as follows:

[0042]

[0043] in:

[0044] c υ is the current growth factor VEGF concentration of the retinal unit grid; the initial value is the steady-state value under normal conditions;

[0045] D is the diffusion coefficient matrix of growth factor VEGF, specifically:

[0046]

[0047] d is a fixed value;

[0048] λ * is the half-life degradation rate, which is a constant value;

[0049] λ is the consumption rate of growth factor VEGF by endothelial cells, which is a constant value;

[0050] η is the production efficiency of growth factor VEGF under the control of oxygen concentration, which is a constant value;

[0051] S(C o ) is the control equation for generating the growth factor concentration under the influence of the current oxygen concentration of the retinal unit grid; if the oxygen concentration of the unit grid is C lim There is no growth factor VEGF generated due to hypoxia. If the oxygen concentration of the unit grid is less than C lim VEGF will be produced due to lack of oxygen, so the equation is as follows:

[0052]

[0053] Optionally, the control equation for the movement direction of the blood vessel tip is as follows:

[0054]

[0055] in:

[0056] is the chemotaxis parameter of growth factor VEGF;

[0057] is the fibrin chemotaxis parameter;

[0058] K is the anisotropic tensor matrix of the extracellular matrix, which is an intrinsic property of the matrix space of each extracellular unit grid (x, y)∈Ω and is independent of time. It is non-uniform and anisotropic in space, which means that the extension direction of the tip cell is not necessarily consistent with the direction of the growth factor concentration gradient or the fibronectin concentration gradient; the control equation is:

[0059]

[0060] k a is the anisotropy parameter, representing the resistance to movement along a certain direction;

[0061] k cond Heterogeneity parameters;

[0062] υ x and y is the component of the unit vector υ that describes the orientation of the extracellular matrix at the unit grid (x, y). is a regular equation, randomly generating υ x and y , and is fixed at all times t; for each unit grid (x, y) ∈ Ω, k cond and k a Also randomly generated.

[0063] Optionally, the control equation of the blood vessel tip movement speed is as follows:

[0064]

[0065] in:

[0066] P is the velocity of the blood vessel tip;

[0067] k p is the rate of vascular tip division, which is related to the concentration of growth factor VEGF c υ Related.

[0068]

[0069] l e is a constant, the length of the endothelial cell;

[0070] t c The cell division cycle is determined by the concentration of growth factor VEGF c υ Determine that the control equation is as follows:

[0071]

[0072] τ is a constant, the cell proliferation time coefficient;

[0073] is the average growth factor (VEGF) concentration; the sum of the growth factor concentrations of all unit grids in the two-dimensional discrete coordinate system at the current time divided by the total number of unit grids, calculated once in each time period t;

[0074] c υ is the current growth factor VEGF concentration;

[0075] c υlim is the threshold concentration of growth factor.

[0076] Optionally, the step S600 includes the following steps:

[0077] According to the velocity of the blood vessel tip in the x and y directions, the vector (x p ,y p ) and the discrete coordinates of the blood vessel tip (x 0 ,y 0 ) Calculate the end point coordinates (x 1 ,y 1 ), where x 1 =x 0 +x p ,y 1 =y 0 +y p ; The movement speed is the first movement speed or the second movement speed;

[0078] The Bresenham algorithm is used to calculate the discrete coordinates (x 0 ,y 0 ) to the end point coordinate of the extension of the blood vessel tip (x 1 ,y 1 ) and the discrete coordinates of the unit grids passed by the straight line are taken as the vascular tip (x 0 ,y 0 ) is a predicted extension distance within this time period t; the predicted extension distance is the first predicted extension distance or the second predicted extension distance.

[0079] Optionally, the conditions for the occurrence of vascular branching when the vascular tip is extended include:

[0080] The blood vessel existence time of the blood vessel tip is greater than a preset blood vessel tip existence time threshold τ a ;

[0081] The unit velocity vector of the blood vessel tip in this time period t The angle change between the angle between the positive direction of the x-axis of the coordinate system and the angle between the unit velocity vector of the blood vessel tip and the positive direction of the x-axis of the coordinate system in the last time period t is greater than the preset angle threshold τ u ;

[0082] When the two conditions of vascular branching are met at the same time, the process of simulating vascular branching at the vascular tip includes:

[0083] Starting from the tip of the blood vessel where the blood vessel branch occurs, the blood vessel is extended in the second movement direction according to the second predicted extension distance and in the first movement direction according to the first predicted extension distance;

[0084] The conditions for vascular fusion to occur when the vascular tip is extended include:

[0085] Among the cell grids included in the blood vessel tip prediction extension distance, there are cell grids whose cell grid attributes are blood vessels and / or blood vessel tips.

[0086] When the vascular fusion conditions are met, the vascular fusion process of simulating the vascular tip includes:

[0087] According to the order from the discrete coordinates of the blood vessel tip to the corresponding end point coordinates on the predicted extension distance, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the predicted extension distance, and extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed;

[0088] The blood vessel tip performing blood vessel extension comprises:

[0089] Each blood vessel tip in the two-dimensional discrete coordinate system is sequentially extended in the current time period t for simulation, wherein the blood vessel tip is extended in the current time period t for simulation as follows:

[0090] Determine whether vascular branching occurs at the vascular tip;

[0091] If the vascular tip does not branch, determine whether the vascular tip has vascular fusion at the first predicted extension distance. If vascular fusion does not occur, set the actual extension distance of the vascular tip to the unit grid to which the vascular tip belongs and all unit grids at the first predicted extension distance, use the current vascular tip existence time of the vascular tip as the initial value of the last unit grid timer at the actual extension distance, reset the vascular tip existence time of the vascular tip to zero and stop timing, change the unit grid attribute of the last unit grid at the actual extension distance to the vascular tip, start the timer to record the vascular tip existence time from the initial value, and change the unit grid attributes of the remaining unit grids at the actual extension distance to the blood vessel, and end the extension simulation of the vascular tip in this t time period. If vascular fusion occurs, the vascular tip is sequentially fused with the blood vessel and the vascular tip at the first predicted extension distance in the order from the discrete coordinates of the vascular tip to the end point coordinates at the first predicted extension distance, and is extended according to the actual extension distance of the vascular tip after the vascular fusion, until the extension simulation of the vascular tip in this t time period is ended.

[0092] If a blood vessel branch occurs at the tip of the blood vessel, the tip of the blood vessel is extended in the second movement direction according to the second predicted extension distance and the tip of the blood vessel is extended in the first movement direction according to the first predicted extension distance; wherein:

[0093] The extending of the blood vessel tip in the first movement direction according to the first predicted extension distance specifically includes:

[0094] determining whether vascular fusion occurs at the tip of the blood vessel at the first predicted extension distance;

[0095] If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the first predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip, and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended;

[0096] If blood vessel fusion occurs, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the first predicted extension distance in the order from the discrete coordinates of the blood vessel tip to the end point coordinates on the first predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after blood vessel fusion until the extension simulation of the blood vessel tip in this time period t is terminated;

[0097] The second movement direction of the blood vessel tip extending according to the second predicted extension distance specifically includes:

[0098] determining whether vascular fusion occurs at the tip of the blood vessel at the second predicted extension distance;

[0099] If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the second predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended;

[0100] If vascular fusion occurs, the vascular tip is fused with the blood vessels and vascular tips on the second predicted extension distance in sequence from the discrete coordinates of the vascular tip to the end point coordinates on the second predicted extension distance, and is extended according to the actual extension distance of the vascular tip after vascular fusion until the extension simulation of the vascular tip in this time period t is terminated.

[0101] Optionally, the vascular tip is fused with the blood vessel and the vascular tip at the predicted extension distance, and the vascular tip is extended according to the actual extension distance of the vascular tip after the vascular fusion, including:

[0102] If the vessel fused with the blood vessel tip is on the predicted extension distance, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel fused with the blood vessel tip belongs on the predicted extension distance, the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the unit grid attributes of all unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip in this time period t is ended;

[0103] If the blood vessel tip that undergoes blood vessel fusion is the blood vessel tip at the predicted extension distance, blood vessel fusion is performed according to the movement direction of the blood vessel tip and the movement direction of the blood vessel tip at the predicted extension distance, specifically:

[0104] If the angle between the moving direction of the blood vessel tip and the moving direction of the blood vessel tip at the corresponding predicted extension distance is greater than or equal to 90 degrees, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel fusion occurs at the predicted extension distance, and the blood vessel tip existence time of the blood vessel tip and the blood vessel tip at the predicted extension distance is cleared and the timing is stopped, and the unit grid attributes of all unit grids at the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip in this time period t is terminated;

[0105] If the angle between the movement direction of the blood vessel tip and the movement direction of the blood vessel tip at the corresponding predicted extension distance is less than 90 degrees, if the blood vessel tip is randomly selected for extension, the blood vessel tip at the predicted extension distance stops extending; when there is no blood vessel or blood vessel tip after the blood vessel tip at the predicted extension distance that merges with the blood vessel tip at the corresponding predicted extension distance, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids included in the predicted extension distance of the blood vessel tip, and the blood vessel tip existence time of the blood vessel tip and the blood vessel tip at the predicted extension distance that merges with it is cleared and the timing is stopped, the unit grid attribute of the last unit grid at the actual extension distance is changed to the blood vessel tip and the timer is started to count from the initial value Record the existence time of the blood vessel tip and change the unit grid attributes of the remaining unit grids on the actual extension distance to blood vessels, and end the extension simulation of the blood vessel tip in this t time period; when there are still blood vessels and / or blood vessel tips after the blood vessel tip on the predicted extension distance that merges with the blood vessel tip on the corresponding predicted extension distance, the blood vessel tip continues to merge with the blood vessel and blood vessel tip on the predicted extension distance, and extends according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this t time period is ended; wherein, the initial value of the timer of the last unit grid on the actual extension distance is the existence time of the blood vessel tip before the blood vessel tip is cleared when no blood vessel branch occurs, and the initial value is zero if a blood vessel branch occurs;

[0106] If the blood vessel tip is randomly selected to stop extending, the blood vessel tip is extended at the predicted extension distance where the blood vessel fusion occurs; the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel tip belongs at the predicted extension distance where the blood vessel fusion occurs on the corresponding predicted extension distance, and the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the attributes of all unit grids on the actual extension distance are changed to blood vessels, ending the extension simulation of the blood vessel tip within this time period t.

[0107] The method for simulating the process of diabetic retinopathy provided by the embodiment of the present invention has the following technical effects:

[0108] (1) The effect of blood glucose levels on diabetic retinopathy is fully considered, and oxygen concentration and blood glucose concentration are incorporated into the simulated vascular extension model. This not only achieves accurate simulation predictions, but also can more comprehensively simulate the occurrence and development of diabetic retinopathy. Understanding how these factors interact to cause retinal vascular abnormalities helps to reveal the root cause of the disease;

[0109] (2) Current treatments for diabetic retinopathy include laser therapy, anti-VEGF injections, etc. By dynamically simulating the diabetic retinopathy process at different stages in this application, it is possible to predict disease progression and the effects of different treatments, optimize treatment strategies at different stages, and improve treatment effects;

[0110] (3) Considering the blood sugar control level and oxygen supply status of individual patients, a personalized simulation model is established for different patients to achieve personalized disease prediction and treatment plan optimization;

[0111] (3) By simulating the regulatory effects of different drugs on oxygen concentration and blood sugar levels, potential new drugs can be screened and their efficacy can be evaluated, thereby accelerating the new drug development process.

[0112] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the written description and the accompanying drawings.

[0113] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0115] Figure 1 A schematic flow chart of a method for simulating a diabetic retinopathy process provided by an embodiment of the present invention;

[0116] Figure 2 Flow chart of step S700 in an embodiment of the present invention. DETAILED DESCRIPTION

[0117] The preferred embodiments of the present invention are described below in conjunction with the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention.

[0118] Oxygen concentration: Studies have shown that retinal hypoxia is a key factor in causing abnormal angiogenesis. However, the specific mechanism of oxygen concentration distribution and its dynamic changes in diabetic retinopathy is still unclear.

[0119] Blood sugar level: High blood sugar is the main feature of diabetes, and it is known that high blood sugar can cause endothelial cell dysfunction and increase oxidative stress. However, the specific mechanism of how blood sugar levels aggravate retinopathy by affecting angiogenesis and oxygen supply still needs further study.

[0120] Comprehensive model construction: This study will be the first to incorporate oxygen concentration and blood glucose levels into the mathematical model of angiogenesis to build a more comprehensive model of diabetic retinopathy. By simulating the retinal angiogenesis process under different oxygen and blood glucose conditions, the interaction mechanism will be revealed.

[0121] In diabetic patients, hyperglycemia is the main feature of diabetes. Hyperglycemia can lead to damage to retinal microvessels and aggravate hypoxia, thus triggering the following series of biological reactions:

[0122] High blood sugar triggers inflammation and cell damage: Long-term high blood sugar can induce abnormal metabolic products in cells, triggering an inflammatory response, leading to retinal tissue cell damage, blood vessel blockage and contraction, and affecting blood supply;

[0123] Exacerbated hypoxia: As high blood sugar affects the normal function of microvessels, retinal tissue hypoxia is further aggravated;

[0124] Release of angiogenic factors: Hypoxia stimulates retinal tissue to release angiogenic factors, such as vascular endothelial growth factor (VEGF), which promotes the growth of new blood vessels;

[0125] Increased angiogenesis: Under the action of angiogenic factors, the retinal tissue begins to abnormally proliferate new blood vessels. These new blood vessels are unstable and prone to rupture and bleeding, seriously threatening vision.

[0126] In summary, in order to simulate the complex phenomena and system behaviors in the process of diabetic retinopathy, help researchers understand and predict the complexity of the process of diabetic retinopathy, analyze how blood glucose levels aggravate retinopathy by affecting angiogenesis and oxygen supply, and the role of oxygen concentration distribution and its dynamic changes in diabetic retinopathy, this application incorporates oxygen concentration and blood glucose concentration into the mathematical model of angiogenesis for the first time to construct a more comprehensive diabetic retinopathy model. By simulating the retinal angiogenesis process under different oxygen and blood glucose conditions, the interaction mechanism is revealed. Figure 1 As shown, the embodiment of the present application provides a method for simulating the process of diabetic retinopathy, the method comprising the following steps:

[0127] Step S100, construct a two-dimensional discrete coordinate system (x, y), and define a vascular index function n(x, y, t) that describes and records the attributes of each retinal unit grid in the two-dimensional discrete coordinate system at each preset time t; the unit grid attributes include blood vessels, blood vessel tips, and environment; it can be seen that the vascular index function n(x, y, t) in the embodiment of the present application is actually a binary value function, with a value of 0 or 1 or 2, indicating the attribute of the unit grid (x, y) at time t, for example, n(x, y, t) is 0, indicating that the unit grid attribute at time t is the environment, and n(x, y, t) is 1, indicating that the unit grid attribute at time t is the environment. When n(x, y, t) is 2, it means that the cell grid attribute is the blood vessel tip at time t. A binary value function is used as the vascular network morphology descriptor, which captures grid information more accurately than the traditional endothelial cell density function. The blood vessel tip is a state before it becomes a blood vessel. The biological processes such as blood vessel tip extension, branching and fusion are used to simulate the generation process of retinal blood vessels. The cell grid attribute is the environment, which simulates the space where blood vessels can grow. The cell grid with the cell grid attribute of the environment may become a cell grid with the attribute of a blood vessel or a blood vessel tip as the blood vessel tip extends, branches and merges.

[0128] Step S200, start simulation for a certain patient, initialize each unit grid attribute, blood sugar concentration, oxygen concentration, fibronectin concentration and growth factor concentration in the two-dimensional discrete coordinate system. In order to be closer to the actual situation of the patient and ensure that the simulation results are closer to the actual development, the blood sugar concentration and oxygen concentration are set according to the patient's examination results, and the initial values ​​are directly obtained from the patient's most recent test results. The initial values ​​of fibronectin concentration and growth factor concentration are steady-state values ​​under normal conditions; the initial values ​​of blood sugar concentration, oxygen concentration, and fibronectin concentration of each unit grid are the same, and the blood sugar concentration, oxygen concentration, and fibronectin concentration of each unit grid change with time. There are differences in blood sugar concentration, oxygen concentration, and fibronectin concentration between each unit grid.

[0129] In an embodiment of the present application, when initializing the attribute of each unit grid in the two-dimensional discrete coordinate system, the unit grid attribute of at least one unit grid is the blood vessel tip, and the attributes of the remaining unit grids are the environment or the blood vessel; the position of the unit grid whose unit grid attribute is the blood vessel tip is not restricted here, but in order to better show the effect of blood glucose concentration on angiogenesis over time, when initializing in this embodiment, the number of unit grids whose unit grid attribute is the blood vessel tip is much smaller than the number of unit grids whose unit grid attribute is the environment, and the number of unit grids whose unit grid attribute is the blood vessel can be zero, or can be less than the number of unit grids whose unit grid attribute is the environment.

[0130] Step S300, at preset time intervals t, update the blood glucose concentration, oxygen concentration, fibronectin concentration, and growth factor concentration of each unit grid according to the current blood glucose concentration of the unit grid, the control equation of the blood glucose concentration, the current oxygen concentration, the control equation of the oxygen concentration, the current unit grid properties, the current fibronectin concentration, the control equation of the fibronectin concentration, the current growth factor concentration, and the growth factor diffusion reaction equation, and save them under the information of the corresponding unit grid.

[0131] In the embodiment of the present application, step S300 includes the following steps:

[0132] The first step is to predict the blood glucose concentration of each unit grid after time t based on the current blood glucose concentration of the retinal unit grid and the control equation of the blood glucose concentration; and save it under the information of the corresponding unit grid.

[0133] The control equation of the retinal blood glucose concentration in the embodiment of the present application is input according to the patient's specific conditions, and is as follows:

[0134]

[0135] in:

[0136] Cg is the current blood glucose concentration of the retinal unit grid; the initial value has been set in step S200;

[0137] C g ( t ) is the blood glucose concentration Cg of the retinal unit grid after time t based on the current blood glucose concentration Cg of the retinal unit grid;

[0138] F ( t ) is the increase in blood glucose concentration at time t after food intake;

[0139] I ( t ) is the increase in insulin concentration after time t, and the initial value of insulin concentration is obtained by patient testing;

[0140] H ( t ) is the increase in glucagon concentration after time t, and the initial value is obtained from the patient's test;

[0141] k 1 is the influence coefficient of insulin concentration on blood glucose concentration, which is a fixed constant;

[0142] k 2 is the influence coefficient of glucagon concentration on blood glucose concentration, which is a fixed constant;

[0143] L ( t) is the blood glucose concentration consumed by the liver effect after time t.

[0144] Step 2: predict the oxygen concentration of each unit grid after time t based on the current oxygen concentration of the retinal unit grid, the blood glucose concentration after time t, and the control equation of the oxygen concentration; and save it under the information of the corresponding unit grid;

[0145] In the embodiment of the present application, the change of oxygen concentration over time is mainly due to the generation, diffusion and consumption of oxygen; in the retina, oxygen mainly comes from retinal blood vessels, and the simulation needs to consider the process of oxygen diffusion from blood vessels to retinal tissue, so the control equation of oxygen concentration in the retina is as follows:

[0146]

[0147] in:

[0148] c o is the current oxygen concentration of the unit grid of the retina, the initial value of which has been set in step S200;

[0149] D o is the diffusion coefficient of oxygen in retinal blood vessels, which is a fixed constant;

[0150] F(c o ) is the oxygen consumption rate equation of the retinal unit grid:

[0151]

[0152] A is the normal oxygen consumption rate of the unit grid, which is a fixed constant;

[0153] B is the constant of oxygen consumption of the unit grid under the influence of blood glucose;

[0154] c olim The lowest oxygen concentration for normal physiological activities of the unit grid is a fixed constant. Unit grids below this constant value are in a hypoxic state.

[0155] G(ρ) is the amount of oxygen transported from the surrounding blood vessels in the unit grid. Vascular transport is the source of oxygen for the retinal unit grid. According to the vascular density at the current time, the calculation equation of G(ρ) is as follows:

[0156]

[0157] α o is the oxygen transport rate constant in blood vessels;

[0158] C is the transport efficiency constant under the influence of blood glucose;

[0159] ρ(x,y,t 0) is the unit grid (x, y) of the retina at the current time t 0 The blood vessel density at the time is specifically the number of cell grids with blood vessels in the current cell grid (x, y) and the nine cell grids adjacent to the current cell grid;

[0160] Step 3: predict the fibronectin concentration of each retinal unit grid after time t according to the current unit grid properties of the retinal unit grid, the current fibronectin concentration and the control equation of the fibronectin concentration, and save it under the information of the corresponding unit grid;

[0161] Fibronectin is attached to the extracellular matrix and is an important substance for the formation of blood vessels. The increase in fibronectin concentration depends on the presence of endothelial cells, and the decrease in fibronectin concentration depends on the uptake of fibronectin by endothelial cells. It can be seen that whether it increases or decreases, it depends on the unit grid whose attribute is a blood vessel, that is, it is related to the value of the blood vessel index function n(x, y, t), that is, t = t 0 When the cell grid property of the cell grid (x, y), the control equation of the fibronectin concentration without the diffusion term in the embodiment of the present application is:

[0162]

[0163] in:

[0164] c f is the current fibronectin concentration of the retinal unit grid; the initial value has been set in step S200;

[0165] ω is the production rate of fibronectin in endothelial cells, which is a fixed constant;

[0166] μ is the uptake rate of fibronectin by endothelial cells and is a fixed constant.

[0167] Step 4: predict the growth factor concentration of each unit grid after time t based on the current oxygen concentration of the retinal unit grid, the growth factor VEGF concentration and the growth factor VEGF diffusion reaction equation, and save the information of the corresponding unit grid;

[0168] The diffusion reaction equation of the growth factor VEGF in the embodiment of the present application is obtained by mass balance including diffusion, accumulation and consumption terms. The situation that the endothelial cells consume the growth factor VEGF (receptor-mediated binding) only occurs when n(x, y, t)=2 and the consumption rate of the endothelial cells for the growth factor VEGF is λ, and is proportional to . When the unit grid attribute is environment (n=0), the growth factor VEGF is generated according to the oxygen concentration; specifically as follows:

[0169]

[0170] in:

[0171] c υ is the current growth factor VEGF concentration of the retinal unit grid; the initial value is the steady-state value under normal conditions;

[0172] D is the diffusion coefficient matrix of growth factor VEGF, specifically:

[0173]

[0174] d is a fixed value;

[0175] λ * is the half-life degradation rate, which is a constant value;

[0176] λ is the consumption rate of growth factor VEGF by endothelial cells, which is a constant value;

[0177] η is the production efficiency of growth factor VEGF under the control of oxygen concentration, which is a constant value;

[0178] S(C o ) is the control equation for generating the growth factor concentration under the influence of the current oxygen concentration of the retinal unit grid; if the oxygen concentration of the unit grid is C lim There is no growth factor VEGF generated due to hypoxia. If the oxygen concentration of the unit grid is less than C lim VEGF will be produced due to lack of oxygen, so the equation is as follows:

[0179]

[0180] Step S400, predicting the movement direction of the blood vessel tip in the two-dimensional discrete coordinate system within this time period t according to the fibronectin concentration of the retinal unit grid after t time, the growth factor concentration after t time and the control equation of the movement direction of the blood vessel tip, and saving it under the information of the corresponding unit grid.

[0181] The dynamics of the vascular tip are influenced by the structure and chemical properties of the extracellular matrix. The direction of motion of the vascular tip u is related to the growth factor concentration c after time t υ and fibronectin concentration c after time t f The gradient of K is proportional to the properties of the extracellular matrix defined by K.

[0182] According to the gradient of growth factor VEGF concentration and fibronectin concentration, the control equation of the direction of blood vessel tip movement is:

[0183]

[0184] in:

[0185] is the chemotaxis parameter of growth factor VEGF;

[0186] is the fibrin chemotaxis parameter;

[0187] K is the anisotropic tensor matrix of the extracellular matrix, which is an inherent property of the unit grid (x, y)∈Ω matrix space outside each cell and is independent of time. It is non-uniform and anisotropic in space, which means that the extension direction of the tip cell (blood vessel tip) is not necessarily consistent with the direction of the growth factor concentration gradient or the fibronectin concentration gradient; the control equation is:

[0188]

[0189] k a is the anisotropy parameter, representing the resistance to movement along a certain direction;

[0190] k cond Heterogeneity parameters;

[0191] υ x and y is the component of the unit vector υ that describes the orientation of the extracellular matrix at the unit grid (x, y). is a regular equation, randomly generating υ x and y , and is fixed at all times t. Similarly, for each position (x, y) ∈ Ω, k cond and k a Also randomly generated.

[0192] Step S500, predict the first movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the first movement direction according to the control equation of the fibronectin concentration, the first movement direction and the movement speed of the blood vessel tip after t time at the blood vessel tip, and predict the second movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the second movement direction according to the control equation of the fibronectin concentration, the second movement direction and the movement speed of the blood vessel tip after t time at the blood vessel tip, and save it under the information of the corresponding unit grid; wherein the first movement direction is the movement direction of the blood vessel tip in this t time period of the blood vessel tip; the second movement direction is the movement direction of the blood vessel tip in the last t time period of the blood vessel tip. In the embodiment of the present application, the control equation of the movement speed of the blood vessel tip is as follows:

[0193]

[0194] in:

[0195] P is the velocity of the blood vessel tip;

[0196] k p is the rate of vascular tip division, which is related to the concentration of growth factor VEGF c υ Related.

[0197]

[0198] l e is a constant, the length of the endothelial cell;

[0199] t c The cell division cycle is determined by the concentration of growth factor VEGF c υ Determine that the control equation is as follows:

[0200]

[0201] τ is a constant, the cell proliferation time coefficient;

[0202] is the average growth factor (VEGF) concentration; the sum of the growth factor concentrations of all unit grids in the two-dimensional discrete coordinate system at the current time divided by the total number of unit grids, calculated once in each time period t;

[0203] c υ is the current growth factor VEGF concentration;

[0204] c υlim is the threshold concentration of growth factor. When the VEGF concentration is higher than this value, the equation starts to work. represents the critical time t c Depends on the current VEGF concentration c υ and the average concentration When VEGF concentration c υ Greater than the threshold c υlim When the critical time t c It will be calculated according to this equation.

[0205] Step S600: According to the discrete coordinates (x 0 ,y 0 ) and the first movement speed to calculate the first predicted extension distance of each blood vessel tip in the first movement direction in the two-dimensional discrete coordinate system, and, according to the discrete coordinates (x 0 ,y 0 ) and the second movement speed are used to calculate the second predicted extension distance of each blood vessel tip in the two-dimensional discrete coordinate system in the second movement direction, and the predicted extension distance is saved in the information of the corresponding unit grid; wherein the predicted extension distance is the distance of the blood vessel tip in the corresponding movement direction and at the corresponding movement speed from the discrete coordinate (x 0 ,y0 ) starts to move and reaches the end point coordinate (x 1 ,y 1 ) passes through all unit grids. In the embodiment of the present application, step S600 specifically includes:

[0206] According to the movement speed of the blood vessel tip (which can be the first movement speed or the second movement speed), the vector (x p ,y p ) and the discrete coordinates of the cell grid corresponding to the blood vessel tip (x 0 ,y 0 ) Calculate the end point coordinates (x 1 ,y 1 ), where x 1 =X 0 +X p ,y 1 =y 0 +y p ;

[0207] The Bresenham algorithm is used to calculate the discrete coordinates (x 0 ,y 0 ) to the end point coordinate of the extension of the blood vessel tip (x 1 ,y 1 ) and use the discrete coordinates of the cell grids that pass through the straight line connected to the blood vessel tip (x 0 ,y 0 ) is a predicted extension distance in the corresponding movement direction (a first predicted extension distance when the movement speed is the first movement speed, and a second predicted extension distance when the movement speed is the second movement speed).

[0208] Step S700: Perform extension simulation on each blood vessel tip in the two-dimensional discrete coordinate system within the current time period t according to the first movement direction of the blood vessel tip, the existence time of the blood vessel tip, the second movement direction of the blood vessel tip, the first movement speed, the second movement speed, the first predicted extension distance, and the second predicted extension distance, and update the unit grid attributes and the existence time of each unit grid in the two-dimensional discrete coordinate system with the extension simulation of each blood vessel tip, and update the simulation time when all blood vessel tips in the two-dimensional discrete coordinate system have completed the extension simulation. The extension simulation includes blood vessel extension, blood vessel branching that occurs when the blood vessel is extended, and blood vessel fusion that occurs when the blood vessel is extended.

[0209] Blood vessel branching refers to the splitting of a blood vessel tip into two new blood vessel tips. This is an important mechanism for the formation of new blood vessels during angiogenesis. When branching occurs, a blood vessel tip (sprout tip) will become two blood vessel tips. In the embodiment of the present application, in order to simulate this process, the embodiment of the present application is defined as follows:

[0210] (I) The conditions for vascular branching are as follows:

[0211] 1. The existence time of the blood vessel tip is greater than a preset blood vessel tip existence time threshold τ a ;

[0212] 2. Unit velocity vector of the blood vessel tip during this time period t The angle change between the angle between the positive (or negative) direction of the x-axis of the coordinate system and the angle between the unit velocity vector of the blood vessel tip (i.e. the original movement direction of the blood vessel tip) and the positive (or negative) direction of the x-axis of the coordinate system in the last time period t is greater than the preset angle threshold τ u .

[0213] (II) When the two conditions of vascular branching are met at the same time, the vascular branching at the vascular tip is simulated as follows:

[0214] A new blood vessel tip is generated from an existing blood vessel tip, that is, the unit grid with the blood vessel tip attribute has two extension attributes of the blood vessel tips, that is, the movement direction of the blood vessel tip corresponding to the unit grid in this time period t is set to the first movement direction and the second movement direction, and starting from the blood vessel tip where the blood vessel branch occurs, it is extended in the original movement direction (second movement direction) according to the second predicted extension distance; and it is extended in the newly predicted movement direction (first movement direction) according to the first predicted extension distance.

[0215] Vascular fusion refers to the meeting and fusion of two vascular tips or one vascular tip with another vascular tip. In the fusion of vascular tip to vascular tip, both vascular tips may disappear, or only one vascular tip may disappear, depending on the specific situation; in the fusion of vascular tip to vascular tip, the vascular tip becomes part of the ring and disappears, no longer extending, that is, the tip no longer exists. In order to simulate this process, the embodiment of the present application is defined as follows:

[0216] (I) The conditions for vascular fusion are as follows:

[0217] Among the unit grids included in the predicted extension distance of the blood vessel tip (including the first predicted extension distance and / or the second predicted extension distance), there are unit grids whose unit grid attributes are blood vessels and / or blood vessel tips.

[0218] (II) When the vascular fusion conditions are met, vascular fusion is simulated at the vascular tip as follows:

[0219] According to the order from the discrete coordinates where the blood vessel tip is located to the coordinates of the end point on the corresponding predicted extension distance (the order of the unit grids and unit grid attributes on the predicted extension distance and the actual extension distance mentioned in this implementation is in this order), the blood vessel tip is sequentially fused with the blood vessel and the blood vessel tip on the predicted extension distance, and extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed, wherein the actual extension distance is all the unit grids (including the unit grid (x0, y0) to which the blood vessel tip belongs) passes from the discrete coordinates (x0, y0) of the blood vessel tip to the end of the extension simulation after the blood vessel branching and / or blood vessel fusion occurs during the movement of the blood vessel tip in the corresponding movement direction and at the corresponding movement speed according to the corresponding predicted extension distance. In the embodiment of the present application, the blood vessel tip is fused with the blood vessel and the blood vessel tip on the predicted extension distance, and extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion includes:

[0220] 1. If the blood vessel fused with the blood vessel tip is on the predicted extension distance (the first predicted extension distance and / or the second predicted extension distance), it means that the blood vessel tip needs to disappear at the unit grid where the blood vessel is located, and the extension of the blood vessel tip will also stop at the unit grid where the blood vessel is located, that is, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel tip belongs on the predicted extension distance (including the unit grid to which the blood vessel belongs), and the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the unit grid attributes of all unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip in this t time period is ended.

[0221] 2. If the blood vessel tip that undergoes blood vessel fusion is the blood vessel tip at the predicted extension distance (the first predicted extension distance and / or the second predicted extension distance), blood vessel fusion is performed according to the movement direction of the blood vessel tip and the movement direction of the blood vessel tip at the predicted extension distance, specifically:

[0222] (1) If the angle between the movement direction of the vascular tip (which may be the first movement direction or the second movement direction) and the movement direction of the vascular tip at the corresponding predicted extension distance (which may be the first predicted extension distance or the second predicted extension distance) is greater than or equal to 90 degrees, it is considered that both tips become inactive. At this time, the vascular tip where the vascular fusion occurs and the vascular tip at its predicted extension distance will disappear, that is, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all unit grids before the unit grid to which the vascular tip belongs where the vascular fusion occurs at the predicted extension distance, and the vascular tip existence time of the vascular tip and the vascular tip at the predicted extension distance is cleared and the timing is stopped, and the unit grid attributes of all unit grids at the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is terminated.

[0223] (2) If the angle between the moving direction of the blood vessel tip and the moving direction of the blood vessel tip at the corresponding predicted extension distance is less than 90 degrees, it is considered that only one of the two blood vessel tips that have fused has become inactive, and one blood vessel tip is randomly selected to continue extending at the corresponding predicted extension distance, and the other blood vessel tip stops extending at the unit grid where the two blood vessel tips have fused, and no longer extends. Since the two tips are at the same position in the model (i.e., the unit grid where the two blood vessel tips have fused), the selection of which tip becomes inactive will not affect the determinism of the model;

[0224] If the blood vessel tip is randomly selected for extension, the blood vessel tip on the predicted extension distance stops extending; when there is no blood vessel or blood vessel tip after the blood vessel tip on the predicted extension distance that merges with the blood vessel tip on the corresponding predicted extension distance, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids included in the predicted extension distance of the blood vessel tip, and the blood vessel tip existence time of the blood vessel tip and the blood vessel tip on the predicted extension distance that merges with it is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the blood vessel tip and the timer is started to record the blood vessel tip existence time from the initial value, and the other unit grids on the actual extension distance are reset to zero. The unit grid attributes of the remaining unit grids are changed to blood vessels; the extension simulation of the blood vessel tip in this t time period is ended; when there are still blood vessels and / or blood vessel tips after the blood vessel tip in the predicted extension distance that is fused with the blood vessel tip in the corresponding predicted extension distance, the blood vessel tip continues to fuse with the blood vessel and blood vessel tip in the predicted extension distance, and extends according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this t time period is ended; wherein, the initial value of the timer of the last unit grid in the actual extension distance is the existence time of the blood vessel tip before the blood vessel tip is cleared when no blood vessel branch occurs, and the initial value is zero if a blood vessel branch occurs;

[0225] If the blood vessel tip is randomly selected to stop extending, the blood vessel tip is extended at the predicted extension distance where the blood vessel fusion occurs; the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel tip belongs at the predicted extension distance where the blood vessel fusion occurs (excluding the unit grid to which the blood vessel tip belongs at the predicted extension distance), and the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the attributes of all unit grids at the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip within this t time period is ended; regardless of whether the blood vessel tip at the predicted extension distance where the blood vessel fusion occurs completes the extension simulation within this t time period, the blood vessel tip extension simulation continues, and it has no effect on it.

[0226] In the embodiment of the present application, in order to simulate the process of blood vessel extension, the definition is as follows:

[0227] Each blood vessel tip in the two-dimensional discrete coordinate system is sequentially extended in the current time period t by simulation, and the blood vessel tip is extended in the current time period t by simulation including:

[0228] Determine whether vascular branching occurs at the vascular tip;

[0229] If the vascular tip does not branch, determine whether the vascular tip has vascular fusion at the first predicted extension distance. If vascular fusion does not occur, set the actual extension distance of the vascular tip to the unit grid to which the vascular tip belongs and all unit grids at the first predicted extension distance, use the current vascular tip existence time of the vascular tip as the initial value of the last unit grid timer at the actual extension distance, reset the vascular tip existence time of the vascular tip to zero and stop timing, change the unit grid attribute of the last unit grid at the actual extension distance to the vascular tip, start the timer to record the vascular tip existence time from the initial value, and change the unit grid attributes of the remaining unit grids at the actual extension distance to the blood vessel, and end the extension simulation of the vascular tip in this t time period. If vascular fusion occurs, the vascular tip is sequentially fused with the blood vessel and the vascular tip at the first predicted extension distance in the order from the discrete coordinates of the vascular tip to the end point coordinates at the first predicted extension distance, and is extended according to the actual extension distance of the vascular tip after the vascular fusion, until the extension simulation of the vascular tip in this t time period is ended.

[0230] If a blood vessel branch occurs at the tip of the blood vessel, the tip of the blood vessel is extended in the second movement direction according to the second predicted extension distance and the tip of the blood vessel is extended in the first movement direction according to the first predicted extension distance; wherein:

[0231] The extending of the blood vessel tip in the first movement direction according to the first predicted extension distance specifically includes:

[0232] determining whether vascular fusion occurs at the tip of the blood vessel at the first predicted extension distance;

[0233] If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the first predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip, and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended;

[0234] If blood vessel fusion occurs, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the first predicted extension distance in the order from the discrete coordinates of the blood vessel tip to the end point coordinates on the first predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after blood vessel fusion until the extension simulation of the blood vessel tip in this time period t is terminated;

[0235] The second movement direction of the blood vessel tip extending according to the second predicted extension distance specifically includes:

[0236] determining whether vascular fusion occurs at the tip of the blood vessel at the second predicted extension distance;

[0237] If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the second predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended;

[0238] If vascular fusion occurs, the vascular tip is sequentially fused with the blood vessels and vascular tips on the second predicted extension distance in the order from the discrete coordinates of the vascular tip to the end point coordinates on the second predicted extension distance, and is extended according to the actual extension distance of the vascular tip after vascular fusion until the extension simulation of the vascular tip in this time period t is terminated.

[0239] like Figure 2 As shown, in this embodiment, step S700 includes the following steps:

[0240] Step S701, judging whether a blood vessel branch occurs in the extension simulation of the blood vessel tip within this time period t according to the first movement direction of the blood vessel tip, the existence time of the blood vessel tip, the second movement direction of the blood vessel tip and the blood vessel branching condition, if a blood vessel branch occurs, execute step S711; otherwise, execute step S721.

[0241] Step S711, the blood vessel tip is extended according to the first predicted extension distance and the second predicted extension distance respectively, and then steps S712 and S713 are executed;

[0242] Step S712, judging whether vascular fusion occurs when the vascular tip is extended according to the first predicted extension distance according to the unit grid attributes of each unit grid on the first predicted extension distance and the vascular fusion condition; if vascular fusion occurs, executing step S714; otherwise, executing step S715;

[0243] Step S713, judging whether vascular fusion occurs when extending according to the second predicted extension distance according to the unit grid attributes of each unit grid on the second predicted extension distance and the vascular fusion condition, if vascular fusion occurs, executing step S716; otherwise, executing step S717;

[0244] Step S714: in the order from the discrete coordinates of the blood vessel tip to the coordinates of the end point on the first predicted extension distance, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the first predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed; then step S702 is executed;

[0245] Step S715, setting the actual extension distance of the blood vessel tip to the unit grid to which the blood vessel tip belongs and all the unit grids on the first predicted extension distance, clearing the blood vessel tip existence time of the blood vessel tip and stopping the timing, changing the unit grid attribute of the last unit grid on the actual extension distance to the blood vessel tip and starting the timer to record the blood vessel tip existence time from zero, and changing the unit grid attributes of the remaining unit grids on the actual extension distance to blood vessels, and ending the extension simulation of the blood vessel tip in this time period t; then executing step S702;

[0246] Step S716, in the order from the discrete coordinates of the blood vessel tip to the coordinates of the end point on the second predicted extension distance, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the second predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed; then step S702 is executed;

[0247] Step S717, setting the actual extension distance of the blood vessel tip to the unit grid to which the blood vessel tip belongs and all the unit grids on the second predicted extension distance, clearing the blood vessel tip existence time of the blood vessel tip and stopping the timing, changing the unit grid attribute of the last unit grid on the actual extension distance to the blood vessel tip and starting the timer to record the blood vessel tip existence time from zero, and changing the unit grid attributes of the remaining unit grids on the actual extension distance to blood vessels, and ending the extension simulation of the blood vessel tip in this time period t; then executing step S702;

[0248] Step S721, determining whether blood vessel fusion occurs at the tip of the blood vessel at the first predicted extension distance, if blood vessel fusion does not occur, executing step S722; otherwise, executing step S723;

[0249] Step S722, in the order from the discrete coordinates of the blood vessel tip to the coordinates of the end point on the first predicted extension distance, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the first predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed; then step S702 is executed;

[0250] Step S723, setting the actual extension distance of the blood vessel tip to the unit grid to which the blood vessel tip belongs and all the unit grids on the first predicted extension distance, and using the current blood vessel tip existence time of the blood vessel tip as the initial value of the last unit grid timer on the actual extension distance, clearing the blood vessel tip existence time of the blood vessel tip and stopping the timing, changing the unit grid attribute of the last unit grid on the actual extension distance to the blood vessel tip and starting the timer to record the blood vessel tip existence time from the initial value, and changing the unit grid attributes of the remaining unit grids on the actual extension distance to blood vessels, ending the extension simulation of the blood vessel tip in this time period t, and then executing step S702;

[0251] Step S702, determine whether all blood vessel tips in the two-dimensional discrete coordinate system have completed the extension simulation within the current time period t. If all have completed, execute step S703; otherwise, extend the next blood vessel tip in the two-dimensional discrete coordinate system within the current time period t and execute step S701.

[0252] Step S703, the simulation time T plus the preset time t, and save the unit grid properties of each unit grid in the two-dimensional discrete coordinate system after time t, and then execute step S704.

[0253] Step S704, determine whether the simulation end condition is met, if the simulation end condition is met, output the simulation result and end, otherwise, enter the next extension simulation of the time period t, and repeat the above process from step S300; wherein the simulation result includes but is not limited to the simulation time T, the initial value of the patient's blood glucose concentration, the initial value of the oxygen concentration, the unit grid attribute of each unit grid in the two-dimensional discrete coordinate system after the last extension in the time period t, and the unit grid attribute of each unit grid in the two-dimensional discrete coordinate system after each extension in the time period t. The simulation end condition can be a preset simulation time threshold, or it can be the simulation time T input by the user according to the needs, or the percentage of the unit grids in the two-dimensional discrete coordinate system that are blood vessels to the total number of unit grids in the two-dimensional discrete coordinate system, and the embodiment of the present application does not limit this.

[0254] A method for simulating the process of diabetic retinopathy provided in an embodiment of the present application fully considers the impact of the blood glucose level of diabetic patients on diabetic retinopathy, and incorporates oxygen concentration and blood glucose concentration into the simulated vascular extension model, which not only simulates the occurrence and development process of diabetic retinopathy more comprehensively and understands how these factors interact to cause retinal vascular abnormalities, which helps to reveal the root cause of the disease, but also achieves accurate simulation prediction of diabetic retinopathy, providing a reliable basis for optimizing and adjusting diabetic retinopathy treatment plans and developing and verifying new drugs; in addition, the embodiment of the present application also considers the blood glucose control level and oxygen supply status of individual patients, establishes personalized simulation models for different patients, and realizes personalized disease prediction and treatment plan optimization.

[0255] Obviously, those skilled in the art can make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims

1. A method for simulating diabetic retinopathy process, characterized in that: The following steps are involved: Step S100, constructing a two-dimensional discrete coordinate system, and defining a blood vessel index function n(x, y, t) that describes and records the unit grid attributes of each unit grid in the two-dimensional discrete coordinate system at each preset time t; wherein the unit grid attributes include blood vessels, blood vessel tips, and environment; Step S200, initializing the blood glucose concentration and oxygen concentration of each unit grid in the two-dimensional discrete coordinate system according to the patient's examination results, and initializing the unit grid properties, fibronectin concentration and growth factor concentration of each unit grid in the two-dimensional discrete coordinate system according to a preset initial state; Step S300, updating the blood glucose concentration, the oxygen concentration, the fibronectin concentration, and the growth factor concentration of each unit grid at every preset time t according to the current blood glucose concentration of the unit grid, the control equation of the blood glucose concentration, the current oxygen concentration, the control equation of the oxygen concentration, the current unit grid properties, the current fibronectin concentration, the control equation of the fibronectin concentration, the current growth factor concentration, and the growth factor diffusion reaction equation; Step S400, predicting the movement direction of each blood vessel tip in the two-dimensional discrete coordinate system within the current time period t according to the fibronectin concentration of the unit grid after t time, the growth factor concentration after t time, and the blood vessel tip movement direction control equation; Step S500, predicting the first movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the first movement direction according to the control equation of the fibronectin concentration at the blood vessel tip after t time, the first movement direction and the movement speed of the blood vessel tip, and predicting the second movement speed of each blood vessel tip in the two-dimensional discrete coordinate system when it moves in the second movement direction according to the control equation of the fibronectin concentration at the blood vessel tip after t time, the second movement direction and the movement speed of the blood vessel tip; wherein the first movement direction is the movement direction of the blood vessel tip in the current t time period of the blood vessel tip; and the second movement direction is the movement direction of the blood vessel tip in the last t time period of the blood vessel tip; Step S600, calculating a first predicted extension distance of each blood vessel tip in a two-dimensional discrete coordinate system in a first movement direction according to the discrete coordinates (x0, y0) of the blood vessel tip and a first movement speed, and calculating a second predicted extension distance of each blood vessel tip in a second movement direction according to the discrete coordinates (x0, y0) of the blood vessel tip and a second movement speed; wherein the first predicted extension distance or the second predicted extension distance is all unit grids through which the blood vessel tip moves in the first movement direction or the second movement direction at the first movement speed or the second movement speed from the unit grid discrete coordinates (x0, y0) to which the blood vessel tip belongs, and reaches the end point coordinates (x1, y1) after a preset time t; Step S700, performing extension simulation on each blood vessel tip in the two-dimensional discrete coordinate system within this time period t according to the first movement direction of the blood vessel tip, the existence time of the blood vessel tip, the second movement direction of the blood vessel tip, the first movement speed, the second movement speed, the first predicted extension distance and the second predicted extension distance, and updating the unit grid attributes and the blood vessel tip existence time of each unit grid in the two-dimensional discrete coordinate system with the extension simulation of each blood vessel tip, and updating the simulation time when all blood vessel tips in the two-dimensional discrete coordinate system complete the extension simulation; wherein the extension simulation includes blood vessel extension, blood vessel branching occurring during blood vessel extension and blood vessel fusion occurring during blood vessel extension.

2. The method according to claim 1, characterized in that The step S300 includes the following steps: Predict the blood glucose concentration of each unit grid after time t according to the current blood glucose concentration of the unit grid and the control equation of the blood glucose concentration; and save it under the information of the corresponding unit grid; Predict the oxygen concentration of each unit grid after time t based on the current oxygen concentration of the retinal unit grid, the blood glucose concentration after time t, and the control equation of the oxygen concentration; and save it under the information of the corresponding unit grid; Predicting the fibronectin concentration of each retinal unit grid after time t according to the current unit grid properties of the retinal unit grid, the current fibronectin concentration and the control equation of the fibronectin concentration, and saving it under the information of the corresponding unit grid; According to the current oxygen concentration, growth factor concentration and growth factor diffusion reaction equation of the unit grid, the growth factor concentration of each unit grid after time t is predicted, and the information of the corresponding unit grid is saved.

3. The method according to claim 1, characterized in that: The governing equation for the oxygen concentration is as follows: in: c o is the current oxygen concentration of the retinal cell grid; D o is the diffusion coefficient of oxygen in retinal blood vessels, which is a fixed constant; F(c o ) is the oxygen consumption rate equation of the retinal unit grid: A is the normal oxygen consumption rate of the unit grid, which is a fixed constant; B is the constant of oxygen consumption of the unit grid under the influence of blood glucose; c olim The lowest oxygen concentration for normal physiological activities of the unit grid is a fixed constant. Unit grids below this constant value are in a hypoxic state. G(ρ) is the amount of oxygen transported from the surrounding blood vessels in the unit grid. Vascular transport is the source of oxygen for the retinal unit grid. According to the vascular density at the current time, the calculation equation of G(ρ) is as follows: α o is the oxygen transport rate constant in blood vessels; C is the transport efficiency constant under the influence of blood glucose; ρ(x,y,t0) is the blood vessel density of the retinal unit grid (x,y) at the current time t0, specifically the number of unit grids with blood vessels in the current unit grid (x,y) and the nine unit grids adjacent to the current unit grid.

4. The method according to claim 1, characterized in that The control equation of the fibronectin concentration is: in: c f is the current fibronectin concentration of the retinal unit grid; the initial value has been set in step S200; ω is the production rate of fibronectin in endothelial cells, which is a fixed constant; μ is the uptake rate of fibronectin by endothelial cells and is a fixed constant.

5. The method according to claim 1, characterized in that The growth factor diffusion reaction equation is as follows: in: c v is the current growth factor VEGF concentration of the retinal unit grid; the initial value is the steady-state value under normal conditions; D is the diffusion coefficient matrix of growth factor VEGF, specifically: d is a fixed value; λ * is the half-life degradation rate, which is a constant value; λ is the consumption rate of growth factor VEGF by endothelial cells, which is a constant value; η is the production efficiency of growth factor VEGF under oxygen concentration control, which is a constant value; S(C o ) is the control equation for generating the growth factor concentration under the influence of the current oxygen concentration of the retinal unit grid; if the oxygen concentration of the unit grid is C lim There is no growth factor VEGF generated due to hypoxia. If the oxygen concentration of the unit grid is less than C lim VEGF will be produced due to lack of oxygen, so the equation is as follows:

6. The method according to claim 1, characterized in that The control equation for the movement direction of the blood vessel tip is as follows: in: is the chemotaxis parameter of growth factor VEGF; is the fibrin chemotaxis parameter; K is the anisotropic tensor matrix of the extracellular matrix, which is an intrinsic property of the matrix space of each extracellular unit grid (x, y)∈Ω and is independent of time. It is non-uniform and anisotropic in space, which means that the extension direction of the tip cell is not necessarily consistent with the direction of the growth factor concentration gradient or the fibronectin concentration gradient; the control equation is: k a is the anisotropy parameter, representing the resistance to movement along a certain direction; k cond Heterogeneity parameters; v x and v y are the components of the unit vector v that describes the orientation of the extracellular matrix at the cell grid (x, y). is a regular equation, randomly generating v x and v y , and is fixed at all times t; for each unit grid (x, y) ∈ Ω, k cond and k a Also randomly generated.

7. The method according to claim 1, characterized in that The control equation of the blood vessel tip movement speed is as follows: in: P is the velocity of the blood vessel tip; k p is the rate of vascular tip division, which is related to the concentration of growth factor VEGF c v Related; l e is a constant, the length of the endothelial cell; t c The cell division cycle is determined by the concentration of growth factor VEGF c v Determine that the control equation is as follows: τ is a constant, the cell proliferation time coefficient; is the average growth factor (VEGF) concentration; the sum of the growth factor concentrations of all unit grids in the two-dimensional discrete coordinate system at the current time divided by the total number of unit grids, calculated once in each time period t; c v is the current growth factor VEGF concentration; c vlim is the threshold concentration of growth factor.

8. The method according to claim 1, characterized in that The step S600 includes the following steps: According to the velocity of the blood vessel tip in the x and y directions, the vector (x p ,y p ) and the discrete coordinates (x0, y0) of the blood vessel tip to calculate the end point coordinates (x1, y1) of the current extension of the unit grid; where x1 = x0 + x p , y1=y0+y p ; The movement speed is the first movement speed or the second movement speed; The Bresenham algorithm is used to calculate the discrete coordinates of the unit grids on the straight line connecting the discrete coordinates (x0, y0) of the blood vessel tip to the end point coordinates (x1, y1) of the current extension of the blood vessel tip, and all the discrete coordinates of the unit grids passed are used as the predicted extension distance of the blood vessel tip (x0, y0) in this time period t; the predicted extension distance is the first predicted extension distance or the second predicted extension distance.

9. The method according to claim 8, characterized in that The conditions for the occurrence of vascular branching when the vascular tip extends include: The blood vessel existence time of the blood vessel tip is greater than a preset blood vessel tip existence time threshold τ a ; The unit velocity vector of the blood vessel tip in this time period t The angle change between the angle between the positive direction of the x-axis of the coordinate system and the angle between the unit velocity vector of the blood vessel tip and the positive direction of the x-axis of the coordinate system in the last time period t is greater than the preset angle threshold τ u ; When the two conditions of vascular branching are met at the same time, the process of simulating vascular branching at the vascular tip includes: Starting from the tip of the blood vessel where the blood vessel branch occurs, the blood vessel is extended in the second movement direction according to the second predicted extension distance and in the first movement direction according to the first predicted extension distance; The conditions for vascular fusion to occur when the vascular tip is extended include: Among the unit grids included in the blood vessel tip prediction extension distance, there are unit grids whose unit grid attributes are blood vessels and / or blood vessel tips; When the vascular fusion conditions are met, the vascular fusion process of simulating the vascular tip includes: According to the order from the discrete coordinates of the blood vessel tip to the corresponding end point coordinates on the predicted extension distance, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the predicted extension distance, and extended according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this time period t is completed; The blood vessel tip performing blood vessel extension comprises: Each blood vessel tip in the two-dimensional discrete coordinate system is sequentially extended in the current time period t for simulation, wherein the blood vessel tip is extended in the current time period t for simulation as follows: Determine whether vascular branching occurs at the vascular tip; If the vascular tip does not branch, determine whether the vascular tip has vascular fusion at the first predicted extension distance. If vascular fusion does not occur, set the actual extension distance of the vascular tip to the unit grid to which the vascular tip belongs and all unit grids at the first predicted extension distance, use the current vascular tip existence time of the vascular tip as the initial value of the last unit grid timer at the actual extension distance, reset the vascular tip existence time of the vascular tip to zero and stop timing, change the unit grid attribute of the last unit grid at the actual extension distance to the vascular tip, start the timer to record the vascular tip existence time from the initial value, and change the unit grid attributes of the remaining unit grids at the actual extension distance to the blood vessel, and end the extension simulation of the vascular tip in this t time period. If vascular fusion occurs, the vascular tip is sequentially fused with the blood vessel and the vascular tip at the first predicted extension distance in the order from the discrete coordinates of the vascular tip to the end point coordinates at the first predicted extension distance, and is extended according to the actual extension distance of the vascular tip after the vascular fusion, until the extension simulation of the vascular tip in this t time period is ended. If a blood vessel branch occurs at the tip of the blood vessel, the tip of the blood vessel is extended in the second movement direction according to the second predicted extension distance and the tip of the blood vessel is extended in the first movement direction according to the first predicted extension distance; wherein: The extending of the blood vessel tip in the first movement direction according to the first predicted extension distance specifically includes: determining whether vascular fusion occurs at the tip of the blood vessel at the first predicted extension distance; If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the first predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip, and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended; If blood vessel fusion occurs, the blood vessel tip is sequentially fused with the blood vessels and blood vessel tips on the first predicted extension distance in the order from the discrete coordinates of the blood vessel tip to the end point coordinates on the first predicted extension distance, and the blood vessel tip is extended according to the actual extension distance of the blood vessel tip after blood vessel fusion until the extension simulation of the blood vessel tip in this time period t is terminated; The second movement direction of the blood vessel tip extending according to the second predicted extension distance specifically includes: determining whether vascular fusion occurs at the tip of the blood vessel at the second predicted extension distance; If the vascular fusion does not occur, the actual extension distance of the vascular tip is set to the unit grid to which the vascular tip belongs and all the unit grids on the second predicted extension distance, and the vascular tip existence time of the vascular tip is cleared and the timing is stopped, the unit grid attribute of the last unit grid on the actual extension distance is changed to the vascular tip and the timer is started to record the vascular tip existence time from zero, and the unit grid attributes of the remaining unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the vascular tip in this time period t is ended; If vascular fusion occurs, the vascular tip is fused with the blood vessels and vascular tips on the second predicted extension distance in sequence from the discrete coordinates of the vascular tip to the end point coordinates on the second predicted extension distance, and is extended according to the actual extension distance of the vascular tip after vascular fusion until the extension simulation of the vascular tip in this time period t is terminated.

10. The method according to claim 9, characterized in that The vascular tip is fused with the blood vessel and the vascular tip at the predicted extension distance, and the vascular tip is extended according to the actual extension distance of the vascular tip after the vascular fusion, including: If the vessel fused with the blood vessel tip is on the predicted extension distance, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel fused with the blood vessel tip belongs on the predicted extension distance, the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the unit grid attributes of all unit grids on the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip in this time period t is ended; If the blood vessel tip that undergoes blood vessel fusion is the blood vessel tip at the predicted extension distance, blood vessel fusion is performed according to the movement direction of the blood vessel tip and the movement direction of the blood vessel tip at the predicted extension distance, specifically: If the angle between the moving direction of the blood vessel tip and the moving direction of the blood vessel tip at the corresponding predicted extension distance is greater than or equal to 90 degrees, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel fusion occurs at the predicted extension distance, and the blood vessel tip existence time of the blood vessel tip and the blood vessel tip at the predicted extension distance is cleared and the timing is stopped, and the unit grid attributes of all unit grids at the actual extension distance are changed to blood vessels, and the extension simulation of the blood vessel tip in this time period t is terminated; If the angle between the movement direction of the blood vessel tip and the movement direction of the blood vessel tip at the corresponding predicted extension distance is less than 90 degrees, if the blood vessel tip is randomly selected for extension, the blood vessel tip at the predicted extension distance stops extending; when there is no blood vessel or blood vessel tip after the blood vessel tip at the predicted extension distance that merges with the blood vessel tip at the corresponding predicted extension distance, the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids included in the predicted extension distance of the blood vessel tip, and the blood vessel tip existence time of the blood vessel tip and the blood vessel tip at the predicted extension distance that merges with it is cleared and the timing is stopped, the unit grid attribute of the last unit grid at the actual extension distance is changed to the blood vessel tip and the timer is started to count from the initial value Record the existence time of the blood vessel tip and change the unit grid attributes of the remaining unit grids on the actual extension distance to blood vessels, and end the extension simulation of the blood vessel tip in this t time period; when there are still blood vessels and / or blood vessel tips after the blood vessel tip on the predicted extension distance that merges with the blood vessel tip on the corresponding predicted extension distance, the blood vessel tip continues to merge with the blood vessel and blood vessel tip on the predicted extension distance, and extends according to the actual extension distance of the blood vessel tip after the blood vessel fusion, until the extension simulation of the blood vessel tip in this t time period is ended; wherein, the initial value of the timer of the last unit grid on the actual extension distance is the existence time of the blood vessel tip before the blood vessel tip is cleared when no blood vessel branch occurs, and the initial value is zero if a blood vessel branch occurs; If the blood vessel tip is randomly selected to stop extending, the blood vessel tip is extended at the predicted extension distance where the blood vessel fusion occurs; the actual extension distance of the blood vessel tip is set to the unit grid to which the blood vessel tip belongs and all unit grids before the unit grid to which the blood vessel tip belongs at the predicted extension distance where the blood vessel fusion occurs on the corresponding predicted extension distance, and the blood vessel tip existence time of the blood vessel tip is cleared and the timing is stopped, and the attributes of all unit grids on the actual extension distance are changed to blood vessels, ending the extension simulation of the blood vessel tip within this time period t.

Images