In the blood flow of the human body, the size of immune cells is generally between 7 and 15 μm (Luo et al., 2011). To simplify the calculations, we use a cell model with a diameter of 10 μm for benchmark simulations. The cell model consists of viscous incompressible cell fluid and linear elastic cell membrane. The cell membrane thickness is 0.4 μm, the elastic modulus is 100 Pa, and the Poisson’s ratio is 0.25 (Aubry et al., 2015). Here, we choose the elastic modulus range of 50–500 Pa to reflect the physiologically possible mechanical response of the cell membrane. Cell fluid and plasma were designated as homogeneous fluids with the same parameters, with a kinematic viscosity of 0.0008 Pa s. The density of immune cells and plasma is 1,000 kg/m3. The phenomenon of cells passing through pores generally occurs in the capillary or venule vascular environment, where the blood flow rate is generally controlled at Re between 0.01 and 0.2 (Haber et al., 2013). Therefore, Re = 0.1 was adopted as the baseline plasma flow field. To investigate the retention effect of microvascular stenosis on cell transport, we selected six sets of parameters with channel widths of 10, 12, 14, 16, 18, and 20 μm to establish a microvessel model containing stenotic pores (Vollmar et al., 2019).

The immune cell transmigration model established in this paper includes a cell droplet and a microvascular model with narrow pores. The initial position of the immune cell is located on the upper side of the pore, and it will passively deform and pass through the pores smaller than its size under the drive of the fluid. Although a 3D model would be more accurate, we chose a 2D model to simplify the calculations. Studies have shown that the model is insensitive to depth in cells entering narrow pores (Leong et al., 2011). At the same time, an axisymmetric model containing half of the computational domain was adopted to analyze the transmigration of cells. The hemisphere model composed of the blue and cyan areas in the figure is the immune cell model, composed of a uniform cytosol wrapped by a linear elastic film. Both cellular fluid and extramembrane plasma are considered homogeneous incompressible Newtonian fluids. This paper regards the plasma flow in microvessels as Poiseuille flow in circular conduits. The velocity of the flow field near the wall is relatively small, and the velocity from the near wall to the microvessel axis rises parabolically and reaches a maximum value. The velocity of the flow field is symmetrical about the microvascular axis. Immune cells are located at the axis of the microvessels and are subjected to a symmetrically distributed fluid shear force. In this paper, COMSOL’s “laminar flow” is used to quantify the fluid flow of cell fluid and plasma. Immune cells are elastically deformed under the drive of blood flow, and “solid mechanics” is used to calculate the deformation of the cell membrane. We used “Fluid-Structure Interaction” to calculate the force-displacement transfer at the fluid-structure interface. The rust-colored area in the figure is the blood vessel wall. This paper ignores its deformation and regards it as a rigid wall.

The CFD model in this paper mainly uses COMSOL multiphysics simulation software for numerical calculation. The governing equation used in the model calculation is the Navier-Stokes (N-S) equation, ρ f [ ∂ U f ∂ t + ( U f . ∇ ) U f ] = ∇ [ μ ∇ U f + μ ( ∇ U f ) T ] + F − ∇ P (1)

The boundary condition is a no-slip boundary, namely the blood flow velocity near the wall is zero. The PARDISO solver is used in the COMSOL solver setting, which uses the stiffness decomposition matrix inversion to solve the N-S equation and uses parallel computing to improve the solution efficiency. The Euler method is used for meshing in COMSOL, and the motion of the fluid-solid boundary is described using an arbitrary Lagrangian-Euler (ALE) method. Distortion of the mesh near the fluid-solid boundary can be solved by remeshing. The COMSOL multiphysics simulation software has meshing capabilities that meet the computational accuracy requirements, and all meshing is completed within the software.

The ALE coordinate system, independent of the Euler and Lagrange coordinate systems, is generally not completely fixed in space, nor is it completely fixed on the material node. And the ALE grid could achieve the appropriate motion.

The Navier-Stokes equation described by ALE is ρ ( v t [ x ] + ( c · ∇ x ) v ) = ∇ x · σ + ρ b (2) The mesh velocity described by ALE can be easily obtained v ^ = v − ∂ x ∂ Χ w (3) where w is the velocity of the solid.

The model is calculated in a two-dimensional environment, and the leukocytes can move along the axial and radial directions of the tube with two degrees of freedom.

The numerical precision of the calculation adopts the default settings of the software. CFD tasks are completed on different computers due to too many cases to be solved.

The so-called Poiseuille flow is the flow of viscous fluid in a circular pipe. When the Reynolds number is less than 2000, the liquid flow in a straight circular pipe with equal cross-section is laminar flow. Poiseuile law Q =  π r 4 Δ p / ( 8 η L ) (4)

It describes the steady flow of an incompressible viscous fluid in a horizontal circular pipe, and the Reynolds number is not large. When the flow form is laminar flow, the flow rate Q and the pressure difference Δp at both ends of the pipe, the pipe radius r ₀ , the pipe length l and the fluid, the relationship of viscosity coefficient η. Poiseuille’s law is an essential law of fluid dynamics. Poiseuille flow is one of the few cases where there is an analytical solution to the N-S equation, the x-axis is chosen to be on the vessel axis, and y is the radial coordinate measured outward from the tube axis. Both the circumferential and radial velocity components are zero, and the velocity component parallel to the tube axis is denoted as u, which depends only on y. At the same time, the pressure is constant in each cross-section. In this way, in the N-S equation expressed in cylindrical coordinates, only one axial equation is left, which simplifies to: μ ( d 2 u d y 2 + 1 y d u d y ) = d p d x (5)

According to the no-slip boundary condition, y = R ,   u = 0 . Available velocity distribution: u ( y ) = − 1 4 μ d p d x ( R 2 − y 2 ) (6)

It can be seen that the velocity from the proximal wall to the microvascular axis increases parabolically.

At y = 0, i.e., on the tube axis (the center of the vessel) u m = R 2 4 μ ( − d p d x ) (7)

The basic unit of the neural network algorithm is the single-layer perceptron model. The structure of the single-layer perceptron model is simple, and its components only have an input layer and a single output layer. The figure below represents a perceptron model with two input neurons and one output neuron node. x ₁ and x ₂ represent the stimulus to the neuron, that is, the input information. w ₁ and w ₂ represent the connection weights between input and output nodes. b represents the excitation threshold of the neuron, that is, the bias. y is the output of the perceptron. The information transfer formula of the single-layer perceptron is expressed as follows: y = ω 1 x 1 + ω 2 x 2 + b (8)

The perceptron works by finding a linearly separable hyperplane separating the loaded data. What’s more, stronger data plane classification capabilities can be obtained by integrating multiple perceptron models. In the original perceptron model, the essence of the inter-layer transfer of the neural network is to perform a linear transformation of the data. This determines that a simple perceptron model can only solve linearly separable problems, but cannot deal with non-linear problems. Based on this, applying a non-linear activation function to the neuron node endows the neural network with the ability to solve linear inseparable problems. The non-linear activation function completes the linear transformation to fulfill output characteristics. As shown in Figure 1, the figure represents the non-linear modification of the linear transformation. The application of non-linear activation functions introduces non-linear factors into the neural network, allowing it to classify data planes using smooth surfaces. This enables the neural network to have superior data classification and regression capabilities.

Assume that the number of nodes in the input, hidden, and output layers are l, m and n. The weights and biases from the input layer to the hidden layer and from the hidden layer to the output layer are in order ω_ij, ω _jk and a _j , b _k , learning rate η. The excitation function is g(x). The detailed steps of BPNN modeling are as follows:

BPNN is a multi-layer feedforward neural network based on an error back propagation algorithm. The principle is to use the common gradient descent algorithm and gradient search technology to minimize the error mean square error between the actual output value and the expected output value of the network in successive update iterations. The model topology of BPNN consists of an input layer, a hidden layer, and an output layer. Figure 2 shows a simple three-layer neural network. The input signal of the neural network acts on the output node through the hidden layer node and generates the output signal. By adjusting the connection weights of input nodes and hidden layer nodes, the connection weights and thresholds between the hidden layer and output nodes during the training, the deviation between the actual output and the expected output value decreases along the gradient direction. After repeated training and parameter adjustment, the network parameters (connection weights and thresholds) that match the minimum error are determined. The network parameters determined by the trained neural network process the input information of the input samples with similar characteristics and output the prediction information. Then the prediction result with the smallest error can be obtained.

Assume that the number of nodes in the input, hidden, and output layers are l, m and n respectively. The weights and biases from the input layer to the hidden layer and from the hidden layer to the output layer are w _ij , w _jk and a _i , b _k . The learning rate is η and the excitation function is g(x). The detailed steps of BPNN modeling are as follows:

(1) Network initialization. First, the weights and thresholds between the input and the hidden layers are initialized. Among them, the common activation function is the sigmoid function. Its function expression is,

Assuming that Y _k is the expected output of the network, the error calculation formula is E = 1 2 ∑ k = 1 n e k 2 (12) e k = Y k − O k (13)

(4) Update of weights and biases. The weights and biases are updated through the back-propagation of the error to make the error function reach the minimum value, where the update formulas of the weights and biases are respectively

Furthermore, we utilize mean squared error (MSE) to measure the accuracy of BPNN predictions. The formula for calculating the MSE is, M S E = S S E l = 1 l ∑ i = 1 n ω i ( y i − y i ∧ ) (16)

The network training of BPNN includes the following processes: (1) Integrate the training data set. The training of BPNN requires large-scale data support, and the training samples need to contain both input data and label data. The accuracy of the data directly affects the prediction accuracy of the neural network. Data normalization is often necessary for the specific training to obtain similar weighting tendencies for different features. Normalization can reduce the computational complexity and speed up the convergence, and also facilitate the subsequent training optimization because the output of the activation function is 0-centered.

(2) Network structure design. The network structure aims to clarify the number of network layers, hidden layer nodes, the activation function between layers, the training and the loss function. As the network structure is determined, all data dimensions are established.

We used the CFD data of immune cell transmigration as the training data of the BPNN model. Considering the limited effect of immune cell stiffness on cell motion and deformation and flow field in the transmigration simulation, it is not used as an input parameter of the neural network. This paper aims to construct a time-dependent prediction model for leukocyte transmigration to establish a model of immune cell movement and flow field. The data we use are all from the data obtained by our CFD simulation, including displacement parameters, cell deformation parameters, flow velocity, fluid lift, and drag parameters. We have fully investigated the hydrodynamic parameters of the environment in which immune cells migrate to ensure the accuracy of the data obtained by CFD. We performed blood flow Re of 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, and 0.2 and AR of 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, a total of 121 CFD simulations. Among them, the combination of Re and AR is (0.1, 0.5), (0.1, 0.55), (0.11, 0.5), (0.12, 0.5), (0.13, 0.5) and (0.14, 0.5). The impaction phenomenon is not included in the training. We will exclude 115 sets of CFD simulation data from these six sets of examples as the basic data set for training.

The 115 groups of Re and AR parameters in the basic data set are used as input data. The transmigration is the target vector to train the time-dependent leukocyte transmigration prediction model based on the BPNN model. In addition, we used 14,221 samples consisting of Re, AR, and cell downflow displacement parameters as a training dataset to train cell motion and flow field prediction models during transmigration. Among them, the target parameters are the cell deformation, the downstream velocity, the fluid lift, and drag forces. In order to test the accuracy of the model, 10 groups of Re and AR were used as (0.2, 0.53), (0.19, 0.57), (0.18, 0.63), (0.17, 0.67), (0.16, 0.73), (0.15, 0.77). New examples of (0.14, 0.83), (0.13, 0.87), (0.12, 0.93) and (0.11, 0.97) are used as test datasets to validate the time-dependent transmigration prediction model. The above parameter combination simulates the migration environment in which most immune cells are located and guarantees the generality of this study. At the same time, four groups of leukocyte transmigration CFD examples with Re and AR of (0.16, 0.73), (0.15, 0.77), (0.14, 0.83) and (0.13, 0.87) were used to test the accuracy of the prediction model of the transmigration. The hydrodynamic simulation still uses the leukocyte transmigration CFD model shown in Figure 3.

Many studies have shown that a 4-layer neural network with only two hidden layers can predict most regular data matrices accurately. The input parameters of the time-dependent leukocyte transmigration prediction model are two vectors, Re and AR, so the number of neurons in the input layer is 2. To avoid the long training time caused by many nodes, we draw on the empirical formula for node selection in the hidden layer, ∑ n = 0 M C I n > h (17) Where M is the number of nodes in the input layer, namely 2. I is the number of hidden layer nodes to be used. h is the number of training samples of the BP model. After calculation, the number of nodes in the hidden layer should not be less than 15. Therefore, we use two-layer hidden layers with 5 and 10 nodes, respectively, to construct the time-dependent cell transmigration prediction model. The activation function from the input layer to the hidden layer and the hidden layer adopts the tansig function, and the activation function from the hidden layer to the output layer adopts the purelin function. In addition, the trainlm function is used as the training function of the BP neural network. Meanwhile, the maximum training times, training target and learning rate are 1,000, 1 × 10^–5 and 0.03, respectively. Figure 4 shows the training performance of the time-dependent leukocyte transmigration prediction model. The model achieves the best verification performance when the number of training times is 38, and its value is 4.524 × 10^–5.

Based on the default partitioning method, the BPNN model randomly assigns 115 samples to the training set, validation set and test set according to the ratio of 60%, 20%, and 20%. Figure 5 shows that the regression accuracies on the training set, validation set, test set, and total dataset reached 0.99959, 0.99907, 0.99874, and 0.99949, respectively. This further proves the prediction robustness of the BPNN model.

To test the prediction accuracy of the leukocyte transmigration time prediction model, we predicted the transmigration time of 10 groups of cases in the test data set different from the training samples and compared them with the CFD results. We use the coefficient of determination as a criterion to analyze the deviation between the BPNN and CFD simulation and further measure the robustness of the mechanical prediction system. The formula for calculating the coefficient of determination is: R 2 = ∑ ( y j ^ − y − ) 2 ∑ ( y j − y − ) 2 (18) Where y j ^ is the predicted value of BPNN and y − is the average value of the real value of the CFD simulation, y _j is the true value of the CFD simulation. The interval range of the coefficient of determination is [0, 1], and the closer to 1, the better the prediction effect.

Figure 6 shows the comparison between the predicted time-dependent of BPNN and the time-dependent CFD for 10 groups of leukocyte transmigration examples from the Re interval of [0.1, 0.2] and the AR interval of [0.5, 1]. The results show that the predicted value in the middle region of the AR interval is in good agreement with the true value and slightly worse at the two ends. The R ² was 0.99846, indicating that the leukocyte transmigration prediction model had excellent prediction accuracy. Based on this, we constructed a BPNN-based time-dependent leukocyte transmigration prediction model a three-dimensional map of the predicted leukocyte transmigration time in the Re interval [0.1, 0.2] and AR interval [0.6, 1] (Figure 7). As shown, the leukocyte transmigration time increased significantly with the decrease of Re. Compared with Re, the reduction of AR has a more prominent effect on prolonging the time-dependent leukocyte transmigration. Under the dual influence of Re and AR, the increase in leukocyte transmigration time was more pronounced and manifested as a sharp upward curl on the map.

The prediction model can predict the time-dependent leukocyte transmigration in specific Re and AR intervals. However, the morphological changes of cells during transmigration and the changes in the surrounding flow field remain unclear. The robustness of the BPNN model aided to construct a leukocyte transmigration prediction model. As mentioned above, 14,221 samples including Re, AR, cell down-flow displacement, cell deformation, down-flow velocity, received fluid lift and drag parameters were used as the training dataset to train the BPNN model. The data matrix containing Re, AR and cell forward displacement is used as input data to establish the corresponding relationship with the target data. The target data vectors are R _max , V _y , F _l and F _d, respectively. Using Eq. 14 to calculate the number of neurons in the hidden layer, we use 50 nodes to fill the hidden layer. We construct a two-layer hidden layer with 20 and 30 nodes, respectively. The selection of activation and training functions is the same as the leukocyte transmigration time-dependent prediction model. At the same time, the determination coefficient of Eq. 11 is still used to measure the prediction accuracy of the BPNN prediction model.

First, we train the neural network with cell deformation parameters as target vectors. Unlike the setting threshold of the time-dependent prediction model of leukocyte transmigration, the better prediction accuracy can be obtained by setting the training target as 10^–3. As shown in Figure 8A, the prediction model achieves the best training performance when the number of training times is 34, and its value is 0.00094417. At the same time, a fitting accuracy of 0.99992 was obtained in training on the overall training set with a sample size of 14,221 (Figure 8B).

Next, we use the cell deformation prediction model to predict the cell deformation during the transmigration for the four new examples mentioned above. Their Re and AR are set to (0.16, 0.73), (0.15, 0.77), (0.14, 0.83), and (0.13, 0.87), respectively. Figure 9A–D compares the predicted cell deformation values and the CFD values of the four groups of calculations, respectively. The figure shows the deformation prediction of the cell forward displacement by the prediction model is consistent with the CFD results. The coefficients of determination, which measure the degree of agreement between the BPNN predictions of the four groups of examples and the real solutions of CFD, are 0.9998, 0.9998, 0.9998, and 0.9999, respectively, showing an excellent prediction effect.

By changing the target vector of the prediction model to the velocity vector of the cell forward movement, we constructed a model for predicting the leukocyte transmigration speed. Due to the increased degree of non-linearity between the input data and the target vector, the threshold for the training target was set to 10^–10. The training performance of the model and the data fitting curve are shown in Figure 10A,B. When the number of training times reaches 64, the BPNN model obtains the best training performance with a value of 9.862 × 10^–11. At the same time, the R value of the training data fitting reached 0.99922, indicating a good mapping relationship between the input data and the output data in the forecast model. Subsequently, we tested the established cell forward velocity prediction model. Using the same 4 datasets, the good predictive ability of the forecasting model is revealed. Figure 11A–D all show that the predicted value of BPNN can well fit the CFD solution, and the coefficients of determination of the closeness between the predicted curve and the real curve are 0.9990, 0.9985, 0.9971 and 0.9979.

Based on the same network skeleton of the leukocyte transmigration prediction model, we replace the trained target vector with the cell fluid lift forces to calculate the accuracy of the prediction model in predicting the cell fluid lift forces. Since the non-linearity of the data further increases, the training objective and minimum gradient threshold are set to 1 × 10^–13 and 1 × 10^–10, respectively, to obtain better prediction accuracy. Figure 12A shows that the prediction model achieved the best training performance when the training order reached 153, and its value was 9.970 × 10^–14. At the same time, the fitting curve of the model achieved an accuracy of 0.99867 (Figure 12B).

Similar to the prediction of cell deformation and downstream motion velocity, we also found that the prediction model of the transmigration is highly accurate in predicting the fluid lift force on the cell. Figure 13A–D of show the prediction curves of the prediction model for the fluid lift forces on cells in the four new cases, respectively. The results show that the predicted curve is consistent with the CFD curve, and the difference between the predicted and actual values at each data point is very small. The coefficients of determination of the predicted curve and the real CFD curve in the four groups of examples are 0.9915, 0.9943, 0.9946 and 0.9965, respectively, showing excellent prediction performance.

Finally, we replaced the target vector of the leukocyte transmigration prediction model with the drag force of the cell to observe the prediction ability of the model to the drag force. After adequate iterations were implemented, it was found that the training model with the training target and the minimum gradient threshold set to 1 × 10^–14 and 1 × 10^–12, respectively, could maintain a good prediction of the target. When the number of training times reaches 716, the model has the best training performance with a value of 9.997 × 10^–15 (Figure 14A). At this point, the fitting parameter of the training data reached 0.99093 (Figure 14B). Based on the same test case, we analyzed the prediction accuracy of the prediction model for the drag force of cells during the transmigration. Figure 15A–D show the drag force prediction curves and CFD curves of different examples, respectively. The prediction accuracies of these examples are 0.9257, 0.9451, 0.9512 and 0.9569, respectively. The entire computation process, including CFD and BPNN, is shown in Figure 16. The prediction models have significant discrepancies in the prediction of the drag force in the initial cell movement stage, the pore entry stage and the pore exit stage. The model maintained a high prediction accuracy for drag force in other periods.

Compared with traditional CFD simulation methods, we found that BPNN greatly shortens the simulation time of the model and maintains high accuracy under the condition of sufficient sample size. Once the network model is trained, BPNN can instantaneously predict the time-dependent transmigration of cells. The computer configuration for this article is (Intel(R) Core(TM) i5-4570 3.2 GHz RAM 8G). Based on the same computing resources, the training time of the BP model is less than 10 min in total. In contrast, the total time required for establishing and training the BP model is comparable to that of a CFD example. This fully reflects the punctuality and efficiency of BPNN forecasting.