Calibration of Circulating Tumor Cell’s Model in Narrow Flow

In this article our objective was to calibrate model of Circulating Tumor Cell (CTC). Different types of cancer produce different types of CTCs. For research purposes, we chose to set up the model according to MCF7 breast cancer cell lines, due to the availability of data from laboratory experiments. First, to obtain working model we used mechanics of our already existing RBC model, taking into consideration the differences between RBC and CTC. Next step was to find values for the elastic parameters of the cell model. We have chosen laboratory experiment where the deformability of breast cancer cell passing through narrow microfluidic channel was examined. The channel has similar dimensions as blood capillaries to mimic the in vivo environment. In order to achieve similar behavior of the cell’s model and the real cell in the experiment we needed to set fluid flow according to experimental data. This was achieved by adjusting the fluid force that is dependent on the volumetric flow rate. Due to the long computational time of the simulation we devised a function between the entry time and the volumetric flow rate. Afterwards we also found dependencies between the changes in elastic parameters and entry time, and we were able to set the elastic parameters so they mimic the behaviour from the laboratory experiment. Further work lies in validating these results.


Intoduction
In recent years, the modeling of biological experiments is increasingly becoming more popular. A good model allows, in some situations, to replace large numbers of biologic experiment with computer simulations, saving time and money. In our research group, Cell-in-Fluid, we model fluid flows into which we place deformable objects. This way, we can for example imitate the blood flow. Nearly 40% to 45% of its volume consists of elastic RBCs. Blood also consist from other objects that allow diagnosis to be made when they're detected.
Metastatic disease is the main cause of morbidity and mortality in cancer patients. [1] In a large percentage of these patients circulating tumor cells (CTC) circulate through the blood stream. Approximately 30% to 40% of patients with solid tumors have micrometastases present in bone marrow, blood, or lymph nodes. However, in the early stages of their spread, they can not be effectively detected by conventional clinical methods. Early diagnosis of micrometastatic disease can lead to better prognosis of patients due to earlier and individualized treatment. One way to detect CTC in blood is to use microfluidic devices. [2] The aim of this article was to propose a CTC calibration procedure using a real biology experiment.
For this purpose, we chose a microfluidic channel with a narrowing, which was used in the article [3] to distinguish between healthy and cancer cells based on their physical properties. e-mail: alzbeta.bohinikova@fri.uniza.sk

Our Simulation Model
Our model consists of two parts. We use the Lattice-Bolzman method for modeling the fluid, and the immersed objects are approximated by the Immersed Boundary Method. Fluid is represented by fictitious particles that move among a fixed grid of discrete points. These fictitious particles then transport the information about the speed and the direction in which they are moving during the simulation. At each discrete point, information about how many fluid particles (how dense the fluid is) and how quickly they pass through that discrete point is stored. [4] Immersed objects are represented only by their boundaries, which are formed by a triangular network of points on the membrane. The movement of the individual points depends on the elastic forces acting on them. We can visualize the each connection of two points as a spring. Consequently, if we apply pressure on this spring, its goal is to return to the original relaxed state.
This elastic behavior is ensured by 5 modulus -stretching (preserving relaxed length of the triangle side), bending (preserving relaxed angle between neighboring triangles), local area conservation (preserving the area of each triangle), global area conservation (preserving the global surface of the cell) and volume (preserving the total cell volume). The forces act at every point of the triangular network, and the behavior of the individual points is always the sum of all the forces that act on it.
All simulations are run in ESPResSo open source software, which includes the Object-in-Fluid software package that allows to model elastic objects in the fluid flow. [5]

Volumetric Flow Rate and Fluid Force
When setting flow in simulations, one of the important steps is to set the correct fluid flow. In [6], we verified the linear dependence between the fluid force value and the volumetric flow rate. What we also have to be aware of is how the fluid force setting works in our simulations. Testing different discretizations of the fluid with the same fluid force gave different values of volumetric flow rates, as you can see in Table 1. The reason is that in our simulations, fluid force adjusts the force acting on each discretization point, and therefore the smaller the value of LBgrid, the more points in the discretization, and that translates to greater force acting overall on the fluid. Thus, the volumetric flow rate is greater at the same level of fluid force. As it is sometimes useful to use of the higher LB grid values, we have decided to verify whether there is a linear dependence between the LBgrid change and the fluid force value when adjusting the same flow in the microfluidic channel.
Our expectations were that when the LB grid is 2 times smoother, the volumetric flow rate will be twice as large with the same fluid force. So, if we want to run simulation with the same fluid flow in a microfluidic device, it is enough to set the correct fluid force at different LB grids. We tested this hypothesis on the set of three simulations. We examined whether when LBgrid = 1 the fluid force is x, then with LBgrid = 2 will the fluid force be 2x and with LBgrid = 4 will the fluid force be 4x if we want to achieve the same flow. Results from the simulations are in the Table 2. The results tell us that the relationship between the transition from LBgrid = 1 to 2 is fully working, when we want to set the same volumetric flow rate at LBgrid = 2, we need to increase the fluid force from LBgrid = 1 two times. When transitioning from LBgrid = 1 to 4, the fluid force needs to change 3.5 times and not 4 times. This was verified by running another set of three simulations. Results are in the Table 3  In article [3], the goal was to find the differences between malignant and non-malignant cells. The microfluidic device created for this purpose consists of a narrow channel and a reservoir at each end. The narrow part of this channel has a length of 150µm and its cross section is a square of size 10 × 10m. The shape of the channel can be seen in the Figure 1. In the experiment, the MCF − 10A cell line (healthy cells used as baseline) and the MCF − 7 line (which represents the cancer cells) were used. The difference between these two cells was measured by three parameters: entry time -the time needed for the cell to fully enter into the narrow part of the channel, the transition time -the time from the moment the cell is completely inside the narrow part, until it completely exits and elongation index. Elongation index is defined as EI = D L /D, where D L is the length of the deformed cell inside the narrow part, and D is the initial cell diameter, before entering the narrow part. The most noticeable difference between cancer and healthy cells was measured by the entry time. Cell sizes used in the article ranged from 15 to 23µm . For this calibration, we chose a cell size of 16µm, since it is the size used in other articles and could be used when verifying our calibration.
For this size of the cell, the elongation index was 1.15µm and the cell average velocity in the channel was 175µm . The entry time was shorter for this size compared to the average, it was 0.15 seconds. These are the values by which we want to adjust our elastic coefficients. We pay the greatest attention to the entry time as we can see that the remaining values did not significantly differentiate between the healthy and cancerous cells.

Simulation set up
In this experiment, the cells passed through a relatively narrow channel, resulting in great stress applied on the cell. Running first sets of simulations in order to set the right fluid force, we noticed a rather odd behavior of the cell. At the site of the channel where the fluid velocity was the greatest, the cell was deformed unnaturally. What was unnatural about its deformation was the fact that its membrane passed through itself, which would not happen in real cells in biological experiments. We can see in Figure  2 what this situation looks like. Therefore, it was necessary to add the interaction between the membrane itself in the simulations. We called this interaction a self-cell interaction, and its role is to control how close the points of a membrane are. If they are too close and there is a risk that the membrane will pass through itself, a repulsive force will begin to work.
So far, the same mechanics have been used to repel two cells in simulations. Where repulsive force is defined for two points (two points previously, each from different cell, now two points on the same cell) as a function of their distance, defined in 1.
where d is the distance of the two points on the membrane, d cut is the threshold at which the repulsive force starts to act, a is the scaling parameter, and n determines how steep the reaction is when two points approach each other. We have to keep in mind that neighbouring points should not Table 5. Values measured in biological experiment [3] a n d cut 0,2 1 6 be repelled when they are at their relaxed distance. The resulting parameters are in Table 5. The first step was the calibration of this repellent reaction to stop the membrane from going trough itself. Subsequently, we set all the necessary fluid coefficients so that the fluid flow was the same as the flow from the biological experiment. The relevant fluid force values are in Table 3. As a starting pount we the simulation where our CTC cell is set with an elastic parameter used for RBC cells. Table 6 shows all the necessary parameters needed to run the simulation. Table 6. Simulation parameters Time step 10 −7 s LBgrid 10 −6 m Fluid force 2, 85 · 10 4 N/m 3 Stretching parameter 8 · 10 −6 N/m Bending parameter 3 · 10 −12 Nm Local area parameter 3 · 10 −6 N/m Global area parameter 9 · 10 −4 N/m Volume parameter 5 · 10 2 N/m 2 Radius of the CTC 8 · 10 −6 m Friction parameter 3, 56 · 10 −9 N sm −1 Fluid density 10 3 kg · m −3 Fluid viscosity 1, 5 · 10 −6 m 2 /s However, this simulation had a relatively long calculation time, so we decided to find faster alternatives. Our goal was to run a set of simulations where we will always change only one of the elastic parameters and then, based on the cell entry time values into the narrow part of the channel, determine what range of parameters to test for our simulation. However, this would have taken long. And our long term aim is to use this procedure to calibrate multiple cell types, so we need to be time-effective, even at the cost of some degree of accuracy.

Fluid Force and Entry Time
Therefore, we have decided to determine the relationship between the length of the entry time of the cell to the narrow part and the size of the fluid force. We've run several simulations to get data. These can be seen in Table 6.
We then fitted these data with an exponential function in the form F(x) = ae bx . Where F(x) denotes the entry time, x is the magnitude of the fluid force, and the parameters a and b determine the shape of the curve.
The appropriate values of parameters a and b were found using the Matlab function lsqcurvefit. [7] This function solves nonlinear curve-fitting problems in least-square sense. The optimisation problem is defined as where in our case xdata is the fluid force, and ydata is the observed entry time for given fluid force values. The results were a = 0.0037 and b = −121, 4321. We can see them also in Figure 3. Thus, at fluid force = 0.0000285, which corresponds to the flow of fluid from the biological experiment, the cell entry time with the elastic coefficients that we had set was 0.003687s. Which is faster than the value of 0.15s that arose from the biological experiment. We expected this result because we have not yet calibrated the elastic cell coefficients.

Calibration of parameters
Next step was the calibration of five elastic parameters. We run sets of simulation where we always took one elastic parameter and changed it to a point the simulation would crash. We left the remaining four parameters fixed. The results were, that the two global parameters, volume and global conservation did not have a large impact on the entry time. So they can be set at kag = 0.9 and kv = 0.5. The results for the other three parameters gave us intervals that worked for the simulations. However these gave us only parameter values that do not crush the simulation. Next step was to run more simulations to get dependencies between these three elastic parameters and entry time. We set the fluid force at 0.1 and tested different values. The entry time in s is in Table 9. The dependencies between elastic parameters and entry time can be seen in Figures 4 and 5.   The best fit for these results was linear dependency between the elastic parameters and entry time. However these results were statistically not significant.

Conclusion
In this article we were able to establish the linear dependency between the fluid force and entry time into the narrow part of the microfluidic channel. Such a results helps us with simulations that are computationally long. Similar approach was done with elastic parameters. We derived working elastic parameters for the model. The future work will lie in fine tuning this approach as well as in validating parameters arriving from this method.