A method for vascular growth applied to a microfluidic chip with a topological structure

Abstract

The invention discloses a blood vessel growth method applied to a microfluidic chip with a topological structure, belonging to the field of blood vessel growth simulation. First, the microchannel topology of the microfluidic chip is set, the growth of blood vessels is simulated in the topological structure, the growth direction and growth speed of blood vessel tip cells are judged, and during the growth process of blood vessels, the constraint conditions of the microchannel wall are introduced to strictly simulate the growth process of blood vessels in the microchannel. The blood vessel regeneration simulation carried out on the microfluidic chip with a topological structure constructed by the invention can adapt to the structures of various microfluidic chips, and different blood vessel networks can be formed under various microfluidic chip structures, which can provide a relatively good reference for the setting and pre-simulation of experiments.

Description

Technical Field

[0001] The present invention relates to the field of vascular growth simulation. Background Art

[0002] Existing vascular growth methods all perform growth simulation in the entire growth matrix. Currently, the relatively advanced phase field theory has been applied to the process of vascular growth. Travasso R, E Corvera Poiré, Castro M, et al. Tumor Angiogenesis and Vascular Patterning: A Mathematical Model[J]. PLoS ONE, 2011, 6(5): e19989. Its model can reflect the change of blood vessel diameter. However, the following problems still exist in the construction of the phase field vascular model:

[0003] (1) The governing equation of the phase field vascular model has a fourth-order differential degree, so it is difficult to stably solve this equation. Vilanova et al. used the IGA method for solving, but the solving process is complex and not easy to implement.

[0004] (2) There are very few studies on successfully constructing the phase field vascular model at the 3D level. The 3D vascular morphology constructed by traditional methods is very messy, and it is impossible to clearly distinguish branches and fusions. Therefore, there is a large difference from the actual 3D capillary network.

[0005] And there is no method for vascular growth in a microfluidic chip in the prior art. Summary of the Invention

[0006] Aiming at the problems in the prior art that the calculation of vascular growth simulation is difficult and there is no corresponding technical solution for growth in a microfluidic chip, the present invention proposes a vascular growth method applied to a microfluidic chip with a topological structure, which is the calculation of the growth process of blood vessels in the microfluidic chip.

[0007] The technical solution of the present invention is a vascular growth method applied to a microfluidic chip with a topological structure, and the method includes:

[0008] Step 1: Construct a microfluidic chip with a topological structure;

[0009] Step 2: Calculate the concentration distribution of VEGF, and calculate the gradient distribution of the VEGF concentration according to the concentration distribution of VEGF; the VEGF is a substance that induces vascular growth, and blood vessels always grow in the direction of higher VEGF concentration;

[0010] Step 3: Set the VEGF concentration distribution and the VEGF concentration gradient distribution obtained in Step 2 on the microfluidic chip constructed in Step 1; Initialize the starting vascular buds of the blood vessels. The tip cells at the top of the vascular buds sense the VEGF concentration to determine the growth direction, and initialize i = 0;

[0011] Step 4: Traverse the old tip cells and determine whether the i-th tip cell meets the activation and maintenance conditions. If it meets, then determine whether this tip cell meets the branching conditions; otherwise, go to Step 8;

[0012] The activation and maintenance conditions are as follows: There is a VEGF concentration gradient around the tip cell sufficient for the tip cell to determine the growth direction. Secondly, there is a sufficient endothelial cell density around the tip cell to facilitate the growth of the tip cell; The VEGF concentration gradient and the endothelial cell density are set according to the actual situation; The branching condition is: Calculate the probability factor ξ of branching. If the calculated probability factor is greater than the set threshold, then the branching condition is met; otherwise, the branching condition is not met;

[0013] If the branching condition is met, set the movement flag number to 2; if not, set the movement flag number to 1. The movement flag number represents the number of movement directions;

[0014] Step 5: Calculate the final movement direction of the tip brain cell according to the VEGF concentration distribution at the current position, and determine whether the first grid point extended along the movement direction belongs to the blood vessel domain. The blood vessel domain represents the set composed of the positions occupied by the blood vessels. If so, go to Step 6; otherwise, determine whether the second grid point extended along the movement direction belongs to the blood vessel domain; If it belongs to the blood vessel domain, set the first grid point along the movement direction as the blood vessel domain and no longer update the position of the first grid point to the new tip cell domain, and go to Step 6; otherwise, determine whether the second grid point extended along the movement direction belongs to the chip channel wall; If it belongs to the chip channel wall, go to Step 6; otherwise, calculate the movement speed of the tip cell. If the speed is greater than the set threshold, move two grid points; if the speed is less than the threshold, move one grid point, and then update the position to the new tip cell domain and go to Step 7;

[0015] The method for calculating the movement speed of the tip cell is as follows:

[0016] Among them, represents the movement speed of the i-th tip cell, c is the VEGF concentration gradient, ▽ is the gradient operation, K is the perturbation matrix, and in the two-dimensional case, K is

[0017]

[0018] Among them, θ is the angular range of deviation from the direction determined by various tropic induction effects during the cell movement;

[0019] Step 6: The tip pus cells remain stationary while updating their positions to the new tip cell domain;

[0020] Step 7: Decrement the movement flag count by 1. Check if the movement flag count is equal to 0. If so, proceed to Step 9; otherwise, go to Step 5.

[0021] Step 8: Check if the participation survival time of the tip cell is less than 0. If so, do not update its position to the new tip pus cell domain and proceed to Step 9. Otherwise, decrement the participation survival time by 1, keep the tip cell stationary, update its position to the new tip cell domain, and then proceed to Step 9.

[0022] Step 9: Check if the i-th tip cell is the last tip cell. If so, update the new tip cell domain to the old tip cell domain, set the new tip cell domain to be empty, and proceed to Step 10. Otherwise, increment i by 1 and then go to Step 4.

[0023] Step 10: Check if the time step count has been reached. If so, the technical cell grows; otherwise, go to Step 3.

[0024] The vascular regeneration simulation performed on the microfluidic chip with a topological structure constructed by the present invention can adapt to the structures of various microfluidic chips, and different vascular networks can be formed under various microfluidic chip structures, which can provide a relatively good reference for the experimental setting and pre-simulation. Description of the Drawings

[0025] Figure 1 Schematic diagram of the branching and fusion behaviors during the vascular growth process;

[0026] Figure 2 Schematic diagram of the movement of tip cells;

[0027] Figure 3 Schematic diagram of the movement of tip cells;

[0028] Figure 4 Schematic diagram of a partial process of vascular growth;

[0029] Figure 5 Simulation data in COMSOL and reconstruction diagram in MATLAB;

[0030] Figure 6 Schematic diagram of the evolution of vascular morphology in the microcolumn model, where q is the time step;

[0031] Figure 7 Simulation data in COMSOL and reconstruction diagram in MATLAB;

[0032] Figure 8 Schematic diagram of the evolution of vascular morphology in the branching model;

[0033] Figure 9 Simulation data in COMSOL and reconstructed graph in MATLAB;

[0034] Figure 10 Evolution diagram of blood vessel morphology in a multi-chamber model. Specific implementation mode

[0035] The existing parent blood vessels will perform angiogenesis behavior under the stimulation of exogenous VEGF, that is, blood vessel sprouts grow from a certain point of the parent blood vessel. The top of the blood vessel sprout is a tip cell with pseudopodia. The pseudopodia of the tip cell can sense the VEGF concentration gradient in the environment and lead the blood vessel sprout to grow in the direction of high VEGF concentration. During the process of the tip cell leading the growth of the blood vessel sprout, the tip cell will undergo branching and fusion phenomena. Branching means that a blood vessel sprout splits into two blood vessel sprouts at the top, and fusion means that when the tip cell encounters an existing blood vessel, it will integrate into the existing blood vessel. For the angiogenesis process in the context of microfluidic chip applications, the microfluidic chip can largely reproduce the actual physiological environment of the blood vessel growth process [1]. At the same time, the microfluidic chip can construct various topological structures to generate a rich VEGF concentration distribution to induce blood vessel growth. The existence of its topological structure will have a certain obstructive or guiding effect on the growth of blood vessels, and thus various blood vessel morphologies can be formed.

[0036] Simulation of blood vessel growth under a microfluidic chip

[0037] The technical problems in this paper mainly focus on two aspects. One is the construction of a microfluidic chip with a topological structure, which is mainly realized in COMSOL Multiphysics. After completing the modeling and calculation in COMSOL, the relevant data is imported into MATLAB. The other is the construction of a blood vessel growth model, which realizes the angiogenesis process under the microfluidic chip in MATLAB by combining the calculation results of COMSOL.

[0038] Construction of a microfluidic chip

[0039] For the construction of the geometric shape of a complex microfluidic chip, this paper completes it in COMSOL and calculates the transient distribution of VEGF through the dilute species transport module. Then, the geometric shape information, VEGF concentration, and VEGF concentration gradient distribution are imported into MATLAB and matrixized, and then combined with the blood vessel growth model for coupled simulation. In particular, for the geometric information of the microfluidic chip, it needs to be processed when imported into MATLAB. It is stipulated that the grid points occupied by the microfluidic chip channels (where the tip cells can move) are 0, and the grid points not occupied by the channels are 1. Through this processing, the restrictive effect of the microfluidic chip on the movement of the tip cells can be reflected.

[0040] Construction of a blood vessel growth model under a microfluidic chip

[0041] The angiogenesis simulation draws on the grid model framework of Chaplain and Sun [2-3], but the model in this paper can achieve the effect of blood vessel morphological growth in a more concise way. The Chaplain model originally included partial differential equations for controlling endothelial cell density and VEGF concentration, and the migration of tip cells was achieved by calculating the VEGF concentration gradient at the position of tip cells at each moment. However, in the microfluidic chip, we assume that the VEGF concentration is constant, so there is no need for control equations. For the control of endothelial cells, we make the following treatment: in the network-based angiogenesis model, the growth of blood vessels is achieved by the movement of tip cells occupying the grid points, and the grid points passed by the tip cells become the blood vessel regions. During the movement of tip cells, during the process of tip cells guiding blood vessel growth, branching behavior will occur, that is, multiple tip cells lead blood vessel growth simultaneously. During the movement of tip cells, if they encounter other blood vessels, fusion behavior will occur. The schematic diagram is as Figure 1 shown.

[0042] In Figure 1 , on the left is the existing mother blood vessel, and two blood vessel buds are distributed on the mother blood vessel. The buds start to grow to the right under the induction of the VEGF source on the right. The growth process is mainly completed by the movement of tip cells distributed at the top of the blood vessel buds and the proliferation of endothelial cells following the tip cells. The growth direction of the blood vessel is mainly determined by the movement direction of the tip cells. As the tip cells continue to move, the blood vessel also continues to elongate. Branching behavior will occur during the movement of tip cells, resulting in the branching phenomenon of blood vessels. Branching behavior will produce more tip cells leading blood vessel growth simultaneously. When the tip cells encounter other blood vessels during movement, tip cell fusion behavior occurs, resulting in the fusion phenomenon of blood vessels.

[0043] To explain this complex process from a mathematical perspective, for the set of all tip cells at each moment, it is represented by the following formula

[0044]

[0045] Ω TEC (t) represents the set of all tip cells at time t, Denote the area (lattice point) occupied by the \(i\)-th tip cell, and \(NTEC(t)\) represents the total number of tip cells at time \(t\). During the movement of tip cells, they are considered to be in an active state, and the maintenance of the active state is jointly determined by the VEGF concentration and gradient in the surrounding environment. The VEGF concentration value ensures that tip cells can receive sufficient information from the extracellular environment, and the concentration gradient value ensures that tip cells can determine their moving directions. The position where tip cells are located should also have a sufficient endothelial cell density. However, in the discrete model, it is assumed that the lattice points occupied by tip cells have sufficient endothelial cell density, so for this condition, it is assumed that the model satisfies everywhere at the positions of tip cells.

[0046] Movement and realization of tip cells

[0047] Movement is the most important feature of tip cells and is a determinant of the blood vessel trajectory during blood vessel growth. For the movement of tip cells, the most important thing is to determine the moving direction of tip cells in each time step. Set the moving speed of the \(i\)-th tip cell as

[0048]

[0049] In formula (2), represents the moving speed of the \(i\)-th tip cell, \(c\) is the VEGF concentration gradient, \(\nabla\) is the gradient operation, and this result is obtained by. \(K\) is the perturbation matrix, and \(K\) in the two-dimensional case is

[0050]

[0051] In the formula, \(\theta\) is the angular range of deviation from the direction determined by various chemotactic inductions during cell movement. This setting can ensure the randomness and diversity of cell movement.

[0052] Due to the limitation of grid lattice points, for \(v\) TEC the direction can only be selected from 8, as Figure 2 shown. Therefore, for the \(v\) TEC of each tip cell, it is necessary to transform \(v\) TEC into the lattice point speed through integer transformation. For the magnitude of \(v\) TEC , it is calculated through the dimensional parameter. In the case where the time step is 0.01, that is, 20 min, the magnitude of \(v\) is 30 μm / h, which is consistent with the actual movement of tip cells. represents the position of the \(i\)-th tip cell at time \(t\) n+1 , represents the position of the \(i\)-th tip cell at time \(t\) n , and the relationship between the two is shown in formula (4).

[0053]

[0054] After the movement of the tip cell is completed, it migrates from the lattice point it originally occupied to the next lattice point, and it is assumed that the tip cell that has moved to the new position has the same endothelial cell density as before. Through the change of position, the tip cell has completed its movement within one time step. In the presence of the microfluidic chip, if the lattice point that the tip cell finally decides to move to is a non-channel lattice point, then the tip cell will remain stationary in this time step; otherwise, the tip cell will move to this lattice point. For the tip cells that do not temporarily meet the activation and maintenance conditions, they are immediately inactivated, but wait for the remaining survival time within a certain number of time steps to cope with the possible re-satisfaction of the activation conditions. In the case where the position where it is located meets the activation and maintenance conditions again. If the activation and maintenance conditions are not met within the specified number of time steps, this tip cell is determined to be inactivated and removed from the tip cell domain.

[0055] Capillary Branching and Realization

[0056] To avoid the rigid setting of the branching angle in many studies, that is, when the tip cell meets the branching conditions, it is set to occupy two adjacent lattice points simultaneously to achieve the realization of the branching effect. However, this can also be improved. Consider the branching as two independent movements of the same tip cell at the same time. That is, for the tip cell to branch, its coordinate information will be recorded in the tip cell domain. At this time, there are two identical coordinates in the tip cell domain, one is the coordinate of the original tip cell, and the other is the same coordinate representing the newly branched tip cell. And during the traversal of the tip cell domain, the independent movement determination of these two tip cells will continue. Due to the existence of the perturbation term during the movement process, there is a high probability of different directions, so the branching structure is formed. For each tip cell lattice point in the tip cell domain, in the branching determination at each moment, it is assumed that they all have the same branching probability ξ, and the probability factor ξ takes the value of 0.1.

[0057] , that is, there is a certain probability of performing the branching behavior, and the branching is as Figure 3 shown.

[0058] Capillary Fusion and Realization

[0059] Since the branching process may cause the tip cell to move inside the blood vessel where it is located, and this situation should not be judged as fusion. The fusion behavior generally occurs during the growth of two different blood vessels. Therefore, before the tip cell moves to the area of the existing blood vessel, it must pass through the tissue area without blood vessel coverage. Therefore, after the tip cell moves, the lattice points in one and two lattice point directions of the moving direction of the tip cell are judged: If the first adjacent lattice point in the moving direction is occupied by the blood vessel area, it means that this tip cell is still inside the mother blood vessel that produced it. Since it cannot move inside the blood vessel, this tip cell remains stationary in this situation; If the first adjacent lattice point in the moving direction is a tissue area and is not occupied by the blood vessel, then judge the second lattice point extending in the moving direction. If this lattice point is a tissue area and is not occupied by the blood vessel, then the movement of the tip cell in this time step is completed and no other actions are continued; If the second lattice point extending out is occupied by the blood vessel, it indicates that there is already an existing blood vessel near this tip cell, and at this time, the fusion phenomenon will occur: Set the position of an adjacent lattice point as the blood vessel area, and remove this tip cell from the tip cell domain. At this time, the fusion behavior is completed. If the blood vessel into which the tip cell fuses is the occupied area of the native blood vessel that produced this tip cell, this behavior becomes the self-fusion phenomenon.

[0060] The overall flowchart is as Figure 4 shown.

[0061] We simulate the blood vessel regeneration process under three types of microfluidic chips with different topologies: one is the microfluidic chip with microcolumn structure, which has a relatively obvious impact on blood vessel growth and branching; the second is the microfluidic chip with branching structure, which can guide the growth direction of blood vessels; the third is the microfluidic chip with various chamber sizes. All numerical values in the three simulations are dimensionless parameter values. Some parameters are shown in Table 1.

[0062] Microcolumn model

[0063] The size of the simulation area is 2×2 (2mm×2mm), and the time step is set to 200 steps (3 days). Among them, the original geometric model of COMSOL and the distribution of VEGF and the geometric model and VEGF distribution after MATLAB reconstruction are as Figure 5 shown. The blood vessel morphology evolution diagram is as Figure 6 shown (sampled every 50 steps).

[0064] MATLAB reconstruction can truthfully reflect the geometric distribution and VEGF distribution in COMSOL, and the tip cells can bypass the circular microcolumns and approach the VEGF source during the blood vessel growth process.

[0065] Branching model

[0066] The size of the simulation area is 2×2 (2mm×2mm), and the time step is set to 240 steps (3.6 days). Among them, the original geometric model and VEGF distribution in COMSOL and the geometric model and VEGF distribution after MATLAB reconstruction are as Figure 7 shown. The diagram of blood vessel morphology evolution is as Figure 8 shown (sampling every 60 steps).

[0067] MATLAB reconstruction can faithfully reflect the geometric distribution and VEGF distribution in COMSOL. The simulation shows that blood vessels can form a network according to the specified shape during the blood vessel growth process.

[0068] Multi-chamber model

[0069] The size of the simulation area is 2×2 (2mm×2mm), and the time step is set to 600 steps (9 days). Among them, the original geometric model and VEGF distribution in COMSOL and the geometric model and VEGF distribution after MATLAB reconstruction are as Figure 9 shown. The diagram of blood vessel morphology evolution is as Figure 10 shown (sampling every 150 steps).

[0070] Table 1 Reaction-diffusion blood vessel growth model parameters

[0071]

Claims

1. A method for vascular growth applied to a microfluidic chip with a topological structure, the method comprises: Step 1: Construct a microfluidic chip with a topological structure; Step 2: Calculate the concentration distribution of VEGF, and calculate the gradient distribution of the VEGF concentration according to the concentration distribution of VEGF; the VEGF is a substance that induces vascular growth, and blood vessels always grow in the direction of higher VEGF concentration; Step 3: Set the concentration distribution of VEGF and the gradient distribution of VEGF concentration obtained in Step 2 on the microfluidic chip constructed in Step 1; initialize the starting vascular buds of the blood vessels, and the tip cells at the top of the vascular buds sense the VEGF concentration to determine the growth direction, and initialize i = 0; Step 4: Traverse the old tip cells, and judge whether the i-th tip cell meets the activation and maintenance conditions. If it meets, then judge whether this tip cell meets the branching conditions, otherwise enter Step 8; The activation and maintenance conditions are: there is a VEGF concentration gradient around the tip cell sufficient for the tip cell to determine the growth direction, and secondly, there is a sufficient endothelial cell density around the tip cell to facilitate the growth of the tip cell; the VEGF concentration gradient and endothelial cell density are set according to the actual situation; the branching condition is: calculate the probability factor ξ of branching, and if the calculated probability factor is greater than the set threshold, it meets the branching conditions, otherwise it does not meet the branching conditions; If the branching conditions are met, the movement flag number is set to 2, and if not, the movement flag number is set to 1. The movement flag number represents the number of movement directions; Step 5: Calculate the final movement direction of the tip brain cell according to the concentration distribution of VEGF at the current position, and judge whether the first lattice point extended along the movement direction belongs to the blood vessel domain. The blood vessel domain represents the set composed of the positions occupied by blood vessels. If so, enter Step 6, otherwise, judge whether the second lattice point extended along the movement direction belongs to the blood vessel domain; if it belongs to the blood vessel domain, set the first lattice point along the movement direction as the blood vessel domain, and no longer update the position of the first lattice point to the new tip cell domain, and enter Step 6, otherwise, judge whether the second lattice point extended along the movement direction belongs to the chip channel wall; If it belongs to the chip channel wall, enter Step 6, otherwise calculate the movement speed of the tip cell. If the speed is greater than the set threshold, move two lattice points, and if the speed is less than the threshold, move one lattice point, and then update the position to the new tip cell domain and enter Step 7; The method for calculating the moving speed of the tip cells is as follows: Among them, represents the migration speed of the i-th tip cell, c is the VEGF concentration gradient, ▽ is the gradient operation, K is the perturbation matrix, and in the two-dimensional case, K is Wherein, θ is the angular range of deviation from the direction determined by various tropic induction effects during the movement of the cell; Step 6: The tip cell remains stationary, and at the same time update its position to the new tip cell domain; Step 7: Decrease the movement flag number by 1, and judge whether the movement flag number is equal to 0. If so, enter Step 9, otherwise enter Step 5; Step 8: Judge whether the participation survival time of the tip cell is less than 0. If so, no longer update its position to the new tip cell domain and enter Step 9, otherwise decrease the participation survival time by 1, the tip cell remains stationary, and update its position to the new tip cell domain and enter Step 9; Step 9: Determine whether the i-th tip cell is the last tip cell. If so, update the new tip cell domain to the old tip cell domain, set the new tip cell domain to be empty, and proceed to Step 10. Otherwise, increment i by 1 and then proceed to Step 4; Step 10: Determine whether the time step has been reached. If so, the technical cell grows. Otherwise, proceed to Step 3.

Images