Modeling cell clusters and their near-wall dynamics in shear flow

The studies that compare the metastatic potential of tumor cell clusters in microcirculation to that of single tumor cells show that the clusters contribute significantly to metastasizing. The metastatic potential is conditioned by the presence of the cancer cells near vessel walls. Detailed understanding of dynamical behavior of clusters near the vessel walls can thus elucidate the process of adhesion. We have developed a biomechanical model of cell clusters capable of simulating both strong and weak adhesion among the cells in the cluster in various spatial configurations. We have validated the model on data from cell separation experiments. The developed model has been used to study near-wall dynamics in shear flow with focus on cluster–wall contact. To quantify the presence of cells near walls, we have evaluated metrics involving time of contact and contact area of clusters tumbling and rolling near the wall. The computational results suggest two trends: First, more elastic clusters and clusters of weakly adhesive cells have decreased cluster–wall contact to the walls than rigid clusters or clusters composed of strongly adhesive cells. Second, more spherical cluster shapes tend to drift away from the walls, thus decreasing the wall contact time.


Introduction
It is estimated that a 10 mL sample of blood from the peripheral circulation of a patient with metastatic cancer typically contains up to 100 individual circulating tumor cells (CTCs) and up to 5 CTC clusters among the approximately 50 billion red blood cells, 80 million white blood cells and 3 billion platelets [4]. These estimates vary according to the cancer type and the treatment stage and are typically measured by isolation and enumeration of cancer cells which contributes to the large width of the ranges.
However, it is known that CTC clusters can be 50 times more metastatic than single CTCs [2]. So even though the clusters are observed much less frequently than single CTCs, it is important to understand their ability to travel through narrow capillaries [6] and their ability to seed metastases, i.e., secondary tumors. The mechanisms are not yet fully understood, but one of the hypotheses is that some of the cancer causing cells are protected inside the cluster from the shear stress imposed by the circulatory flow.
It is known that shear stress has an adverse effect on CTCs. The work [24] explored the effects of shear stress on breast cancer cells with different metastatic abilities, cancer cells of ovarian, lung and leukemic origin and found that shear stress of 60dynes/cm 2 (achievable during intensive exercise) kills significant number of CTCs within few hours and prolonged shear stress effectively reduces the viability of highly metastatic and drug resistant breast cancer cells.
In order to cause metastases, the CTC cluster has to come into contact with a blood vessel wall for a sufficiently long time at sufficiently low velocities.
In vessel bifurcations, due to the splitting of the flow a significant number of cells come into contact to the tip of the Y-junction. At the normal conditions, the highest shear rate is right at the tip of Y-junction [13]. Therefore, at this tip, there is a large wall shear rate that may reduce the chance of cells to adhere. However, further down the stream, the wall shear decreases and the cells that are still within short distance from or even touching the wall may touch the wall again and then adhere. It is therefore reasonable to study the behavior of cells and cell clusters under the assumption that their initial position is near the wall or that they are touching the wall.
At the same time, the lift forces are proportional to the shear rate and cause that particles deviate away from the walls [1,12] unless very specific conditions are met [22]. Therefore, it is of interest to carefully study the relatively short but nevertheless critical period of time that the particles spend close to the walls.
In this work, we develop two different cluster models (one with fixed geometry, one with capability to rearrange its cells), consider the effects of elasticity and compare more spherical and less spherical morphologies under various shear rates.
Using these two models, we computationally analyze the contact time and contact areas of CTC clusters in flow at different shear rates. We must emphasize, however, that the adhesion contact under real adhesion may be different due to additional bio phenomena such as activation of specific adhesive proteins. Nevertheless, the size of contact area and the duration of the contact give a sound basis for further investigations. Recent investigations [3] suggest that larger clusters or microemboli with intermediate nonsphericities have the most increased cluster-wall contact. However, this was only done at high shear rates ( 150-1600 s −1 ) and with rigid models with fixed geometry, i.e., individual cells could neither deform nor slide over one another and thus change cluster conformation.
We investigate the migration of clusters away from the contact wall by placing them near the wall. Therefore, the computed contact area depends on initial position of the cluster. Putting the cluster away from the wall is not leading to experience any contacts.

Cell cluster models
Single cell model The models build upon an existing 3D membrane model for single elastic cells immersed in fluid [15] implemented in module PyOIF [16] in an open source scientific computational package ESPResSo [25].
The single cell model relies on a triangular mesh of points that are bonded by five elastic interactions. This mesh models the membrane of the cell. The stretching, bending and local area interactions, which are parameterized by elastic coefficients k s , k b and k al , respectively, model the local elastic properties of a cell membrane, while the global volume and global area interactions, parameterized by k v and k ag , ensure that the total surface and total volume of the cell does not deviate by more than few percent from the values in relaxed state as the cell undergoes deformations caused by flow.
When there are more cells in the simulation, they interact with each other using pointwise non-bonded potentials, e.g., soft-sphere potential, that repulse particles belonging to the Reprinted from [3] with permissions meshes of different cells. Their collective behavior ensures that the cell membranes do not overlap. The same mechanism is applied to cell-wall interactions.
Since we are considering only laminar flows, the fluid calculations are performed using the D3Q19 lattice-Boltzmann (LB) method, which uses a static cubical grid. Each nonboundary grid point contains 19 populations of fictitious particles that propagate and collide with particles from neighboring grid points. The macroscopic fluid properties, such as velocity and density, are then obtained from these populations. For more details about LB method, see [18].
The cell is coupled to the fluid via a bidirectional dissipative immersed boundary method [11]. At each timestep, the elastic forces, fluid forces, cell-cell interactions forces and cell-wall interactions forces are summarized at each cell membrane node and used to propagate the node using the Newton's second law of motion.
The LB method is used in the whole computational domain, inside and outside the cell. However, the inner viscosity of real cells is higher than the outer fluid; for example, for red blood cells, the viscosity contrast between inner and outer fluid is 3 to 1. In our model, the cell contains a suspension of dissipative inner particles with radii 0.5 μm in order to represent the viscosity contrast between the inner and outer cell environments. The effective viscosity of the inner particle suspension achieved via DPD interactions, per our previous work [9], is ∼ 4.5 mPas, with the outer fluid viscosity being 1.5 mPas and density 1000 kg/m 3 .
The DPD particles inside the cell supplement the LBM fluid inside the cell to reach higher viscosity of the fluid. Cluster geometry We based the cluster geometries on the microscopic images of clusters, e.g., Fig. 1.
We consider four basic cluster shapes: doublet (two cells), tri-cluster (three cells forming a rod), L3-cluster (3 cells in an L formation) and tetra-cluster (4 cells placed at the vertices of a regular tetrahedron), shown in Fig. 2. We consider the Fig. 2 Cluster geometries bi, L3, tetra and tri as represented by the two models. Model I in the first row, Model II in the second row doublet (bi-cluster) as the basic shape for which we calibrate the adhesion. The tetra-cluster is relatively spherical, while tri-cluster and L3-cluster represent the less spherical geometries. In all geometries we use spherical cells with radius 7 μm, and contact area radius approximately 2.7 μm (slightly less for the tetra-cluster due to the contact areas being close together). The cell radius and contact area radius were estimated from images of biological cells, as shown in Fig. 1. The basic triangulation of one sphere has 1002 nodes.
We consider two different cluster models. Model I is suitable for modeling strong adhesion and has the meshes of individual cells joined into one cluster mesh. In a simulation, the cells remain connected and while the cluster may deform (depending on the elastic properties of its membrane), it cannot change the relative positions of individual cells. Both models allow also other cluster geometries (e.g., L4 as generalization of L3 or complex clusters with more cells) that are not shown in this work. Model II links the individual cells using a repulsive-attractive potential. In this model, the cells can slide over one another or under large stress even separate from the cluster, so this model is suitable for clusters with weak adhesion. Mechanical properties of cells In both models we consider two different levels of membrane elasticity represented by elastic coefficients for stretching k s , bending k b and conservation of local area k al , summarized in Table 6. The values were selected in such a way that the soft cell clusters are slightly stiffer than red blood cells (for which the model was calibrated in [15]) and stiff cell clusters are almost rigid.

Model I: clusters with strong adhesion
A simple cluster model that does not allow change in morphology and can be used when the cell-cell adhesion is known to be strong, relies on modeling the cells as one unified object. This can be achieved by defining the shape of the cluster and then preparing the corresponding triangulation mesh.
To prepare the triangulation mesh of the whole cluster, we start with multiple identical triangulations of a sphere, each with radius r s . Each triangulation includes a list of nodes coordinates [x, y, z] (their IDs being their order in this list starting at 0) and a list of triangles represented by triplets of integer IDs of nodes. From these lists, we remove nodes and triangles on the neighboring caps and join the neighboring meshes, see Fig. 3. The joining algorithm for two neighboring spheres involves the following steps (the process is repeated for clusters containing more cells): • identify one point of contact on each sphere (not necessarily a mesh node) • select the contact area radius r c • on both spheres remove caps, i.e., remove all mesh nodes closer to their respective points of contact than • identify mesh nodes on the rings around the holes • for each point on the ring, find the closest neighbor that belongs to the ring of the neighboring mesh and replace this pair of nodes with one node which lies halfway between them • change IDs of the nodes in triangles in the second mesh to reflect the fact that they were replaced by new nodes • shift IDs of all triangles in the second mesh to remove the gaps in numbering left by removing the nodes • join the nodes files and the triangles files

Model II: clusters with weak adhesion
Some CTC clusters have the ability to significantly modify their morphology. The cell-cell adhesion appears to be flexible to such extent that when the microfluidic device geometry dictates it, a spherically shaped cluster can transform into a single chain of cells [5].
In order to model this type of behavior, we devised a more flexible connection of the cells, compared to Model I. Instead Joining process for creating a bi-cluster from two single cells: 1. Remove the caps from both spheres, 2. align the two meshes, 3. link the corresponding vertices from the two meshes and replace them with link midpoints, and 4. re-number the IDs of the nodes in the second mesh. The algorithm is repeated to add more cells for more complex clusters of creating fixed bonds between neighboring mesh nodes, this approach takes advantage of the non-bonded interactions. This is the same type of interaction as the one used in our model for repulsing cells from the channel wall or to prevent overlap of single cells in flow. The mechanics of this type of interaction is based on a potential that is defined between two types of particles, i.e., particles of one cell are of one type and the particles of the second cell have another type. Depending on the potential, the interaction can be repulsive, attractive or combined. Lennard-Jones potential For clusters we need an interaction that is both attractive, to hold the cells together, and repulsive, to prevent the overlap of the membranes. To model this type of behavior in coarse-grained simulations, Lennard-Jones potential is often used [23].
The potential is defined as: where r is the distance between the interacting particles, r cut is the cutoff distance and L J scales the strength of the interaction. The parameter σ determines the distance r min = σ 1 6 , a point of minimum of the potential, where the repulsion changes to attraction.
Before applying the adhesion forces, it is important to adjust the spherical mesh to create a flat adhesion surface, as shown in Fig. 4. The first two steps are identical to the steps described for the strong adhesion. Afterward, for strong adhesion, the selected points have to be removed, whereas here they have to be projected onto the adhesion plane. For any type of cluster, the following process is applied to each pair of cells sharing an adhesion surface: • identify the center point of contact area on both neighboring meshes (the center of contact is the midpoint between the two cell centers; center point on each surface is then defined as the mesh point closest to the center of contact) • identify the points of contact (copulas) on both neighboring meshes (the threshold distance from center point is determined from the desired contact radius) • define the adhesion plane for each cell with normal vector defined by the centers of the two cells • project the points on the copulas onto the adhesion plane • adjust the cell positions so that the adhesion surfaces are at r min distance Since the cells are elastic, they are left to relax after the initial geometry adjustment to achieve the steady state between the elastic and the adhesion forces.
Note that it is possible to achieve a symmetrical adhesion surface in some types of clusters. One simple example is the doublet, where one cell is the image and the other is its mirror reflection. This is especially helpful when creating the strong adhesion. When the cells are symmetrical, the mesh nodes on the ring of each cell are perfectly aligned, and it is much easier to join them into one point. This advantage was used for creating the bi-, L3-and tri-cluster. However, this mirroring is not possible for more complex cluster structures, such as the tetra-cluster, that have multiple contact areas per each cells. Adhesion strength The adhesion strength is determined by the size of the contact area and the LJ potential parameters. In the clusters used in this study we set up each individual adhesion between two cells to have the same strength. To determine the strength of a given doublet of cells, we have devised two simulation experiments where the doublet is separated by force, as shown in Fig. 5. In the first experiment, a pulling force is applied to selected mesh points on the outside copulas of the cells, thus modeling a micropipette. (The results were also confirmed by applying the pulling force to all membrane points. In this case, the cells separated at almost the same pulling force, only the shape transitions were different, as expected.) In the second experiment, the cells were separated by flow. More details on both can be found in [10].
Using these separation experiments, the adhesion parameters can be adjusted to achieve the desired adhesion strength.
Considering the cell-cell adhesion measurements of other types of cells found in the literature, 2 − 12n N between human embryonic kidney cells in [14] and 2−5n N for mesoderm and endoderm cells in [19], we set up the cluster which separates when pulling force of 1.6n N is applied, and holds for lower forces. This can be seen in Fig. 7, where the number of contact points can be used as a parameter determining whether the cluster separated (number of contact points is 0) or not. Note that as the applied force increases, the cells separate faster. As the applied force decreases, after the initial adjustment, the contact areas remain constant. Since we were working with two sets of elastic parameters for the cell membrane (soft and stiff), we calibrated the LJ parameters for both sets. The parameters can be found in Appendix A Table 1. We studied the effects of the elasticity on adhesion in [10].

Cell wall contact
The process of cancer cells extravasation through the endothelium is mediated by highly specific receptor-ligand interac-   tions that allow flowing cells to be brought into contact with and firmly adhere to activated endothelial surfaces before penetrating the subendothelial matrix [8]. During this process, a faster cell velocity for example means that adhesion molecules will be brought into contact more frequently, but as the velocity increases the contact time is decreased. Thus, the optimal binding efficiency depends on the contact duration of CTC with endothelial cells, and on the surface area of each cell that is in contact with endothelial cells. The latter two characteristics of a cell-endothelium encounter can be combined to produce a metric called the time integral of contact area that is an indicator for the probability of an extravasation event [20]. To indicate, evaluate and quantify the cluster-wall contact, similarly to our preliminary investigation of cluster models in bifurcations [17], we use the following metrics: t T -total contact time over a selected time period, A contact -variable contact area, A max -maximum contact area, and TICAtime integral of contact area over a selected time period. The variable contact area allows for observation of differences in contacts across the models and their parameters.
The contact area is calculated as follows. We first calculate n c , the number of mesh nodes that are closer to the wall than the range of the cell-wall interaction. Knowing the total surface of the cluster S (in Model I) and total number of mesh points n, we calculate the fraction of the total surface represented by the contact nodes The rotation of the cluster can be inferred from the minimal z-coordinate z min of the cluster. Unlike for sphere-like objects where the minimal z-coordinate does not represent the rotation well, for rodlike objects, e.g., bi-cluster in Fig.  8, z min oscillates. The period of its oscillation is half of the object's rotation period.

Effect of shear rate
The wall shear rate is influenced by the vessel diameter and differs significantly through the microcirculation. In large vessels such as carotid or coronary artery with diameters 2 to 8mm the wall shear rate achieves values between 100 and 500s −1 . In mid-sized mesenteric arteries with diameter 410 − 660μm the wall shear rate achieves values of 170 − 330s −1 [21]. In smaller vessels, e.g., in submucosal veins with diameter of 120μm the shear rate is up to 200s −1 and in capillaries with diameter of 8μm it increases up to 1000s −1 [7]. The overview of physiological shear rate in differentsized vessels is summarized in Table 2.
The reason for considering different cluster models is shown in Fig. 8 inspired by [3]. As mentioned before, in this figure we can see the cluster's rotation represented by z min . If we take two shear rates (400s −1 and 800s −1 ) and rescale time (i.e., taking every other data point from the slower shear rate), we see that for stiff Model I, strong adhesion cluster, Fig. 8 (a), the behavior is quite consistent. Both the z-coordinate of the centroid, z c , and minimal z-coordinate, z min behave in a similar way. This phenomenon was even more pronounced in [3] where they considered completely rigid cells.
If we do the same for a soft Model I cluster, Fig. 8 (b), we see that the correspondence is not so precise and that the soft membrane slows down the overall movement of the cluster. Therefore, using Model I clusters is only appropriate when we know that the adhesion bonds of the modeled cells are very strong.
In cases of low shear rates, allowing sliding of the cells over one another (Model II) alters the flow pattern only for soft clusters, see Fig. 9. In the case of stiff clusters (solid and dashed blue lines in (a) and (b)), the difference between rotations demonstrated by z min is negligible. Both clusters (Model I and Model II) undergo during 0.1s approximately 3.5 periods of z min oscillations. On the other hand, for soft clusters (solid and dashed red lines in (a) and (b)), we can see significant differences. Model I soft cluster (a) undergoes slightly less than three periods while Model II soft cluster (b) undergoes more than 3.5 periods of z min oscillations.
In cases of high shear rates, the cell sliding is evident even for stiff clusters; for more details, see Fig. 16 in Appendix B.
When we consider the contact area over time under shear rates 400s −1 , as shown in Figure 10, which result in cluster velocity approximately 5mm s −1 , we see that the clusters, with the exception of stiff Model I tri-cluster, tend to drift away from the wall. For higher shear rates, this effect is even more pronounced.
In contrast, in simulations with lower shear rates, 100s −1 , simulated for 0.8s, Fig. 11, we see that clusters other than the tetra-cluster come in contact with the wall more. Under this shear rate, the velocity of the cluster is approximately 1.25mm s −1 . The motion of the tri-cluster rod is the most regular, while the asymmetric motion of the L3-cluster causes irregular wall contacts due to its non-symmetrical geometry.
The previous observations can be summed up by stating that clusters made of more elastic cells tend to drift away from the wall and even weak-adhesion cluster composed of stiff cells exhibit this behavior. This suggests a general conclusion that non-rigid clusters-and here non-rigid means clusters made of either weak adhesive cells or soft elastic cells or both-tend to drift away from the walls at high shear rates.
We do not include the number of contacts as a separate metric, because over the simulated time, we only have the opportunity to see a few contacts and the resulting number would be skewed by an arbitrary time frame cutoff. However, comparing the two models of tri-cluster, Fig. 11 left, we see that the Model II (multiple cells) results in slightly higher frequency of rotation and thus more contacts than Model I.

Effect of elasticity
In Fig. 12, we see the contact areas of clusters with varying elasticity with all other parameters held constant. While for Model I the elasticity does not have an impact on the frequency of contacts, there are significant differences for Model II clusters of different elasticity. While the stiff cluster continues to come into contact with the surface regularly,      Fig. 15 Contact area of the same stiff bi-cluster starting in four different initial orientations: x, y, z and general as depicted in Figure 13 in flow with shear rate 100s −1 . a Model I b Model II the soft cluster stops after the second contact. An intermediate elasticity results in more contacts than the soft cluster, but nevertheless, these also cease and the cluster travels away from the wall L3-and tetra-clusters travel away from the wall after few initial contacts in both models that we investigate. Table 3 shows the contact metrics for clusters in shear rate 100s −1 excluding the first initial contact. For soft cells we see that Model II has shorter contact time and smaller contact areas than Model I. For stiff cells the duration of contacts is similar between the two models, while the maximum contact area is slightly smaller in Model II than in Model I. This indicates that as a first approximation, when one wants to model fairly stiff clusters, both models can be used with similar results. However, for soft clusters, the interplay between the membrane elasticity and the adhesion of cells causes differences. And thus an appropriate model should be selected based on the strength of adhesion of cells in the cluster.
Another thing to note in Table 3 are the differences among the cluster geometries. For more spherical shapes (bi-cluster and tetra-cluster) we observe shorter contacts with smaller contact areas compared to less spherical shapes (tri-cluster and L3-cluster). This is true in both models and indicates that the cluster geometry has a significant impact on the cell-wall contacts.
In Appendix B, Fig. 21, we show the heatmaps for the contact areas in the two models for all cluster types. In the heatmaps we see that the soft clusters have larger contact areas and that the model slightly influences the placement of the contact area, but generally speaking the primary determinant of the contact area is the cluster geometrical shape. For better understanding of the placement of the contact areas on the cluster surface, we also provide a video as part of the Supplementary Material.

Effect of initial orientation
Another factor that impacts the cell-wall contacts is the initial orientation of the cluster.
We considered four different orientations as shown in Fig.  13, where in x and y initial positions all cells are in contact with the wall, thus these clusters have larger initial contact area than the z initial orientation. To represent the general orientation, we selected an initial rotation vector [π/4, π/6, 2π/5], which means that the cluster was rotated by π/4 about the x axis, by π/6 about the y axis and by 2π/5 about the z axis.
In Fig. 14, we see the minimal z-coordinate of clusters that started at the four different orientations. The x and z orientations result in the same rolling motion with the x reaching the first peak later due to the fact that the cluster initially lies flat and needs to rotate to position similar to the initial z position. The period and rotational frequency are the same for these two initial positions and almost the same for the two models (0.175s, 5.71s −1 , respectively, in Model I and 0.176s and 5.68s −1 in Model II).
The general rotation exhibits a similar behavior, but stays closer to the wall, up to 3μm (compared to almost 6μm for x and z rotations) with slightly longer periods, (0.182s in Model I, 0.188s in Model II), and thus also slightly lower rotational frequency (5.49s −1 in Model I and 5.32s −1 in Model II).
Note that the y orientation is qualitatively different, with the two cells exhibiting axle-like rolling motion along the wall, which is smoother in Model II.
The difference between the two models can also be seen in Fig. 15, where the y, z and general orientations experience more contacts in Model I than in Model II.

Conclusion
In this work, we have introduced two types of models that can be used to simulate behavior of cell clusters in various flow conditions. Both models are built on the top of our model of elastic cell and use adhesion among the cells within the cluster. Model I connects the cells with bonds that cannot be broken. This approach further explores the ideas from [3], where rigid cells were merged to create one-mesh rigid cluster. We studied how the behavior of more elastic cluster differs from the rigid one. This type of cluster can be used mainly to answer questions about clusters with strong cell adhesion and stable geometry.
However, when one is interested in more complex behavior, more flexibility needs to be added. In line with that, we have introduced Model II which allows the sliding of cells over each other and even separation of cell from one another, which is a behavior that can be observed in vitro, if the geometry of the microfluidic device dictates it.
We have described the construction of the cluster geometries for both models. To model the weak adhesion in Model II, LJ potential was used. By altering the parameters of the potential, the adhesion strength for this model can be adjusted. To do so, we have used two types of separation experiments, one mimicking the micropipette separation and the other simulating elongation flow, where the cluster is separated only by the fluid forces.
With five parameters modulating the elasticity of individual cells within the cluster and with the possibility to modify the three adhesion parameters for each pair of cells in contact, the model is highly adjustable for various cluster types. We have focused on a selection of cluster geometries, specifically bi-, tri-, L3-and tetra-clusters.
Since adhesion is linked with the physical presence of cells near vessel walls, the contact time and the size of the contact area of the CTCs and their clusters with the walls are the key observables defining the cluster-wall contact.
We computed several metrics combining both observables to study the near-wall dynamics.
By varying the shear rate, we found that non-rigid clusters (made of either weak adhesive cells or soft elastic cells or both) tend to drift away from the walls in higher shear rates. This conclusion might be a bit surprising, since the common notion is that very elastic cells stick to the vessel walls. In that case, however, an additional adhesion mechanism is involved based on creation protein bonds between the cell membrane and the vessel wall [8]. In our work, we focus on purely hydrodynamical effects. The reason for such behavior is caused by the ability of more elastic (or less adhesive) cells in a cluster to adapt to the surrounding fluid. This might be caused by minimization of rotation inertia via regrouping to more sphere-like shape.
In lower shear rates of 100s −1 , the soft and weak adhesive clusters migrate away from the walls slower as the lift forces are proportional to the shear rate. Here, the refined metrics has been used to quantify the cluster-wall contact. In general, however, we can conclude that increasing the elasticity decreases the cluster-wall contact by both reducing the total number of contacts and decreasing the contact area of individual contacts of cluster with wall.
Our final observation is linked to the initial orientation of clusters. In case of a vessel bifurcation, a cluster may enter one of the branches in different orientations. We showed that this also influences the consequent cluster-wall contact.
the development of personalized medicine of selected malignant tumor diseases and its impact on life quality," ITMS code: 313011V446, cofinanced by resources of European Regional Development Fund.

Conflict of interest
The authors declare that they have no conflict of interest. On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Simulation parameters
In the following tables, we provide values of simulation parameters. Table 4 lists time and space discretization param-   eters together with fluid physical properties. In Table 5, we provide surface and volume information for the cluster types we used. In Table 6, we list simulation parameters for the cluster membrane. In Table 7, we show parameters for soft sphere, hat and DPD interactions for inner particles, and in Table 8, we provide the parameters of the repulsive soft-sphere interaction between the cell membrane and the channel walls.

B Additional figures
To provide a visual view of the cluster properties, we show snapshots from the simulations. In Fig. 16, we show how cells slide over each other when the cluster moves-this is relevant to Model II only. The initial contact area is colored red, and each cluster is depicted after 0.2s of flow in shear flow with high shear rate. We can clearly see how the contact area shifts. In Figs. 17-20, we provide snapshots of the flow of the four types clusters near wall.
In Fig. 21, we visualize the physical parts of the cluster that effectively touch the channel wall during the rotating movement. The intensity of the coloring represents the time that individual mesh points were in a contact with the channel wall. Fig. 16 Model II, stiff bi-, tri-, L3-and tetra-clusters in flow with high shear rate 1200s −1 shown after almost 0.2s of flow. The initial cell-cell contact areas are colored red. Due to the flow, the cells change their relative positions and while they stay connected into clusters, the contact areas shift     Heatmaps for soft and stiff clusters using both models. The heatmaps include the large initial contact. Note that for better visibility, the color does not represent the same ranges in the different clusters, but does represent the same range for individual cluster types, i.e., the darkest color represents 47ms of cell-wall contact for bi-clusters (first column), 118.4ms for tri-clusters (second column), 36.2ms for L3 clusters (third column) and 19.8ms for tetra-clusters (fourth column) over the 0.8s simulation. We see that the soft clusters have larger contact areas (most notably see, e.g., tetra-cluster in either model) and that the model slightly influences the placement of the contact area (see, e.g., Model I soft bi-cluster vs Model II soft bi-cluster), but generally speaking the primary determinant of the contact area is the cluster geometrical shape