Flow in filopodia: re-organization and the representation of biological entities as computational objects

Abstract

Actin is the most abundant protein in eukaryotic cells. They form filamentous polymers that are organized in different ways within the cell to perform various functions. For instance, prominent parallel bundles of F-actins mediate the formation and dynamics of filopodia that are long, finger-like protrusions of cell membrane occurring in certain cells, like growing neurons. Understanding actin organization dynamics and its regulation is a crucial problem for biologists that cannot be solved exclusively by biological methods, requiring the support of mathematical and computational modelling. In this work, grounded on a previous hypothesis of ours about the cytosol flow within filopodia, we address several modelling challenges posed by the growth of filopodia in neurons. We use alternative stochastic models and particle-centered numerical methods for transport and elongations, as well as an innovative object-oriented modelling-strategy to represent chemical transformations, polymerization, and their regulation. These modelling strategies allowed for simulating elongations 20 times longer than the typical ranges attained by commonly used filopodia diffusion models, and show that our hypothesis is feasible, acting as a proof-of-concept about the importance of considering organization as a key element in biological explanations.

1 Introduction

Biological entities and their organization stem out of physico-chemical elements and interactions. However, biological phenomena reflect perceptions from an ever-changing reality and cannot be explained solely in terms of physical and chemical dynamics (Harold 2005). Biological reality consists of a hierarchy of intertwined (molecular) processes and biological entities that are generally perceived as space–time (4-D) things performing (biological) functions. Functions are only perceived when these entities and processes interact with other elements in their surroundings. On the top of that, cycles and process-circularity are ubiquitous in life phenomena (Kritz et al. 2010). To tame biological complexity, we need to represent its 4-D things as identifiable units and disentangle their inherent organizational status and hierarchical regulation from descriptions of underlying physico-chemical phenomena, giving substance to them and concomitantly isolating the rules of both types of interaction, better understanding emergent aspects and properties. Modelling, particularly computational and mathematical modelling, is central in this quest. However, the elusiveness of addressing circularly dependent changes and several levels of organization at once, as well as, addressing domains of mathematical functions as consequents rather than antecedents pose enormous challenges to modelling biological phenomena.

The cytoskeleton provides a good example of this imperative. Virtually all cell behavior and intra-cellular functions depend on the cytoskeleton, which constantly grows, shrinks, bends and blends in a myriad of ways to support cell activities (Alberts et al. 2002). Its response is often triggered by conditions at large and foreign to the cell. Defective cytoskeletons have multiple consequences and can cause many diseases, such as tumor-cell migration and metastasis, or neural disorders (Birbach 2008). The analysis of cytoskeletal dynamics requires thinking at the interface of biology, biochemistry, biomechanics, and modelling.

At the leading edge of cell membrane, cytoskeleton may be organized into tightly cross-linked parallel bundles of actin polymers (F-actin), which form membrane long finger-like protrusions called filopodia, with the purpose of probing the environment. These prominent bundles consist fundamentally of a 1-D structure whose geometrical characteristics are ruled by molecular processes grounded on physics and chemistry, that essentially mediate their formation, elongation, and spatial dynamics. Hence, they are amenable to be translated into mathematical-computational models before addressing more complicated phenomena like the protrusions’ architecture or how the organization of filaments and their bundles are ruled by (poorly known) biological laws acting at a higher organizational level.

Explaining flow of matter within filopodia is an important biological quest and long-standing scientific challenge. Current models of filopodial dynamics are based on diffusion as the actin transport mode (Atilgan et al. 2006; Baines 2010) and explain only short elongations that are far from representing the observed biological behavior. In a previous work Moura et al. (2016), we addressed diffusion as a G-actin transport mechanism and its ability to provide enough actins to the polymerization process at the tip of the filopodia for various filament lengths, revealing its inadequacy. We did it through simple calculations based on standard models and without using hydrodynamic indicators, because the fluid conditions within filopodia frontally violate the isotropy and homogeneity assumptions necessary for the meaningfulness of this kind of numbers. As a consequence, we advanced a breakthrough hypothesis to illuminate flow of matter in filopodia.

We proposed a model for filopodial fluid dynamics contemplating inward and outward flows intermediated by a fluid re-organization phase induced by the polymerization that occurs at the tip of the filopodia. Thus, the cytosol-actin mixture could flow towards the filopodia tip where the observed polymerization process would re-organize it into a ‘cytosol trapped within an actin polymer bundle’ conformation that brings matter out of the filopodia protrusion. In this previous article we also embraced the view that mathematical and computational modelling should be part of biologists’ intellectual arsenal to dissect and investigate life phenomena, side by side with direct experimental observation in labs.

The mathematical-computational model presented in the sequel is a first realization of ideas introduced in or supporting the hypothesis above. It is a working model whose first aim is to verify the feasibility of the hypothesis proposed in Moura et al. (2016), see Sect. 2.2, Fig. 3, and, second, to illustrate the possibility of building computational models grounded on very basic processes and components that are easily adaptable and reusable from the programming standpoint, allowing for experimentally adding or subtracting biological features (Sect. 3). It aims at understanding preferably to predicting and to maintain a window towards the investigation of open biological questions such as “How is the flow at the tip of the filopodia?”, “How is polymerization of each filament regulated?”, “How does the filopodia bend and grow into different directions?”, “How is the membrane sustained?”, and so on. This test also delineates a minimum set of biological components needed to test the hypothesis, justifying the choices.

Our model for filopodial flow is centered on the re-organization of the actin-cytosol mixture at the tip of the filopodia. It represents only the most basic and experimentally observable processes, events, and objects that underly the dynamics of filopodia-like cell-protrusions, be them biological or physico-chemical. It is important though to remark that these simplifications do not impair or diminish the confirmation of the explanatory potential of our hypothesis for the filopodia flow, provided by the results reported below. The model parameters, already discussed and justified in our 2016 paper (Moura et al. 2016), are reviewed along Sect. 4 or whenever used. Processes underlying life phenomena may be roughly divided into those related to motion and those related to changes, transformations, and entailments of organizations. Organizations in the living interact and are transformed or aggregated but not created. The modules of our model concerned with motion are described in Sect.2 and those with organization in Sect. 3. Section 4 presents some experiments performed on this virtual platform. Conclusions are in Sect. 5.

2 Fluid flow in filopodia: domain and motion aspects

The whole filopodia dynamics occurs in a changing domain whose deformation depends on the behavior of several phenomena inside it. In this first approach, the domain dynamics is highly simplified by supposing its geometry to be almost cylindrical and its deformations to occur along one single direction. Moreover, the domain is divided into compartments and elongate one compartment at a time (Sect. 2.1). All compartments are cylinders except the one at the tip of the filopodia. In this first implementation, the filopodia is supposed to be radially symmetric and only its 2D cross section along the grow-axis is represented (see Sect. 2.2, Fig. 4). The length of the filopodium changes discretely, depending on the particles’ behavior and their presence at the tip-compartment. When a compartment is added, pushing the tip-compartment further, suspended material occupies the enlarged volume following their dynamics but unobstructed by other material in the direction of the tip (see Sect. 2.2, Fig. 5). A similar procedure was employed by Peskin et al. (1993) and Erban et al. (2014).

The cytosol substrate (colloidal water) is not observable and therefore not explicitly represented. Its displacement is inferred from the behavior of the suspended molecules, which are supposed to undergo a Brownian motion with a preferred direction along the filopodia axis. The movement of suspended molecules is kinematically computed from given velocities. Each particle velocity is a function of cytosol diffusion among other factors (Sects. 2.2 and 2.4). In the outgoing flow, there is no suspended molecules in the cytosol trapped in the filament-bundle and, thus, there is no diffusion. The cytosol movement appears explicitly only when the filopodia elongates through the addition of a new domain-compartment, driven by the accumulation of G-actin molecules in the last (tip) domain compartment (Sect. 2.1), and in the outward flow, trapped inside the filament bundle (Sect. 2.3).

2.1 The filopodia flow domain

The domain is two-dimensional and simulates the interior of a filopodium, which appears entirely in the first quadrant (Fig. 1). Here, the radius of a filopodium, $R=100$ nm, is displayed on the y-axis. The x-axis displays the length of a filopodium, which varies with time due to elongation of the microfilament bundle. Fluid enters the filopodium in the region $x = 0$ and $y \in \left[ {0,r_{ - } } \right) \cup \left( {r_{ + } ,R} \right]$, with a constant concentration $C_{0}$. The region $\left[ {0,L_{c} } \right] \times \left[ {r_{ - } ,r_{ + } } \right]$ represents the microfilament bundle, whereby the reorganized fluid outflows from the filopodium. In practice, while the first mentioned region acts as a source of inflowing cytoplasm, the filament bundle is treated as a sink with borders $y = r_{ - }$ and $y = r_{ + }$, which prevents the entry of volume between them.

Fig. 1

(Available at www.biorender.com)

A geometric representation of the modeled domain. Image created by Thiago F. Leal with BioRender

We define that a filopodium has a length $L(t) = L_{C}$. Despite its time dependency, it is implicitly associated with the number of filaments N and interfilament spacing s, and their impact on polymerization, and not only with the elapsed time. As the polymerization occurs, $L_{C}$ increases and eventually displaces the membrane to open enough space for the polymerization to continue. This process is emulated by adapting the idea of a compartmentalized domain and the filament-membrane interaction described elsewhere (Erban et al. 2014; Lan and Papoian 2008; Peskin et al. 1993).

To do so, we define the filopodium length as $L_{C} = n_{C} \cdot h_{C}$, where $n_{C} \in {\mathbb{Z}}_{ + }$ represents the actual number of compartments, and the length of each compartment is $h_{C} = 27$ nm. Compartments will only be opened if all N filaments in the bundle reach the membrane, represented by the right boundary of the most advanced compartment. Polymerization will occur evenly in the bundle allowing that filaments may extent up to same size, and then it can produce membrane deformation. Figure 2 illustrates the compartmentalized model proposed here. It simulates a limiting factor for growth, a fact already predicted (Blanchoin et al. 2014; Erban et al. 2014; Zhuravlev et al. 2010).

Fig. 2

Filopodial compartmentalized model. Filopodium grows from $x = 0$ and has its length represented by a collection of compartments disposed over x-axis. Once every filament in the bundle reaches the right boundary of most advanced compartment, one more is opened so that the filopodium may elongate. Image created by Thiago F. Leal with BioRender

We define $\delta = 2.7$ nm as the discrete fixed step on which the filament tip changes its position, that is, the contribution of a polymerized monomer in filament length. But, to the filopodial tip move forward, all filaments have to reach the same size. Then, N G-actins must be equally polymerized in the N filaments for the tip to move forward. Therefore, in ten steps the tip walks a compartment length $h_{C}$, which means writing $h_{C} = 10\delta$. Then, defining $n_{p}$ as a dynamic count of polymerized monomers, it is possible to compute the distance $l_{p}$ among F-actin bundle tip (red star in Fig. 2) and lower boundary of last compartment, which is:

$$ l_{p} = \frac{{n_{p} \cdot \delta }}{2N}. $$

(1)

The bundle tip’s position changes dynamically while the filopodium elongates and polymerization will always occur at some point $x_{p}$ inside the last compartment. This means that $x_{p}$ is separated from filopodial base a distance of $(n_{C} - 1)$ compartments plus $l_{p}$.

So, we can write:

$$ x_{p} = (n_{C} - 1)h_{C} + l_{p} , $$

(2)

and then

$$ x_{p} = 10 \cdot (n_{C} - 1)\delta + \frac{{n_{p} \cdot \delta }}{2N}. $$

(3)

2.2 Inflow: cytosol (circa) laminar flow and advective transport of G-actins

F-actin in filopodial bundles undergoes constant “treadmilling”, i.e., backflow of the entire F-actin bundle driven by their polymerization at the distal tip, while their concomitant disassembly occurs at the filopodial base. The quantity of polymerization versus disassembly at any point in time will determine the flow rate and the elongation behavior of the organelle. This requires uninterrupted delivery of large amounts of monomeric actin (G-actins) as building blocks to the tip of F-actin bundle. A key challenge to understand filopodial dynamics is the high rate of actin monomers needed through the entire length of these slender and long structures. This poses the intriguing question as to what transport mechanisms may be at play that can deliver G-actin under these conditions in sufficient amounts. Answering this question harbors fundamental explanations for filopodial behavior.

In previous work (Moura et al. 2016), we show that diffusion does not deliver enough G-actin to sustain the backflow. Diffusion seems to be a proper first model for G-actin motion in short protrusions (about or less than 0.5 µm) (Lan and Papoian 2008; Mogilner and Rubinstein 2005), but needs to be combined with other phenomena to explain lengths observed in vivo. As F-actin constant backflow extracts volume from the very tip generating a negative pressure, this triggers a compensatory inflow of cytoplasm. Therefore, the advective effect of the cytoplasm flow is an additional mechanism, as an alternative to complement diffusion. In addition, we propose that polymerization rearrange the fluid into two phases: the inflow cytosol with molecules suspended outflowing as filaments with cytosol trapped between them, as shown in Fig. 3.

Fig. 3

Two-phase model. In the left, the filopodial inflowing suspended molecules rearrange into filaments as long as polymerization occurs and the F-actin bundle retrograde flow consumes the reorganized volume at filaments plus-ends, inducing a compensatory inflow of cytoplasm due to a drop of pressure in the filopodial tip, which generates a volume flux within filopodia. In the right, the reorganization scheme of polymerization inducing inflow of matter into filopodium. In phase A, the cytoplasm, essentially a mixture of cytosol and suspended molecules (G-actins and ABPs), inflows into filopodia towards the tip, where polymerization induces a reorganization of this fluid. In that stage, G-actins are captured by filaments plus-ends through polymerization, which backflows to the body cell, also consuming a volume of cytosol trapped within bundle filaments. It promotes a reorganization of the cytoplasmic fluid, which consists of the model phase B. The reorganization of fluid and its consequent removal of volume in the filopodial tip in phase B induces a compensatory inflow of cytoplasm, back to the phase A. Image by M.V. Kritz

Then, since polymerization is a chemically reactive process, beyond the affinities of G-actins and ABP, the scenario suggests a reaction–diffusion-advection system to be modeled. In this context, molecules may move and react, simultaneously or not, depending on their motions. Hence, a stochastic approach seems feasible and potentially more promising to describe a chaotic system, with a large amount of interactions and mechanisms occurring almost simultaneously. Further, the outflow volume is straightly related to the number of filaments in the bundle N, the spacing s between them, and filopodial length L itself.

2.3 Outflow: the backflow of F-actin filaments’ bundle

We consider that the membrane is sufficiently rigid to support filopodial elongation and consequent drop of pressure in the very tip. Although the development of a diffusion–advection model in this context should consider the inflow and outflow volumes in the filopodium and their respective velocities. The inflowing fluid is essentially cytosol containing a number of G-actins or profilin-actin complexes. The outflow actin filament bundle is hexagonally packed (Fig. 4B) (Grazi 1997; Jasnin et al. 2013) and the filaments are cross-linked by specific proteins, such as fascin (Harker et al. 2019).

Fig. 4

Hydration shell around actin filaments. A Schematic representation of a single microfilament with its hydration coat, which is composed of water molecules attracted to the polar actins. B Schematic example of a cross-linked bundle, where every microfilament and cross-linkers are surrounded by shells of hydration, filling a fraction of inter-filament spacing and increasing the filaments diameters. Notice that this may leave some spaces unfilled. Image created by Thiago F. Leal with BioRender

Cross-linkers are located every 25–60 actin subunits and bring mechanical cohesion and stiffness to the bundle (Blanchoin et al. 2014; Grazi 1997; Mogilner and Rubinstein 2005), also suppressing any lateral movement of these actin filaments. Several authors have analyzed the relation between the length and/or binding strength of cross-linkers and the resistance of bundles to buckling (Mogilner and Rubinstein 2005; Zhuravlev and Papoian 2011). Then, based on those works, there is a theoretical foundation to argue that inter-filament spacing is narrow, fairly constant, and varies in the range of few nanometers.

The most accepted values for inter-filament spacing were mentioned to range from 7 to 12 nm (Islam 2004). However, filaments can be expected to be hydrated; actin molecules are polar and their charges can attract water molecules (Kabir et al. 2003). This is a chemical property of both G-actin and F-actin. In the case of actin filaments, hydration is a cellular self-regulatory process to compensate high osmotic stress. Once intracellular osmotic control is reestablished, filaments' hydration layers become less thick (Matsuo et al. 2016; Schwienbacher et al. 1995). This directly impacts on the structure of the filament bundles. The hydration layer increases filaments diameters, thus increasing the outflowing volume and reducing the unoccupied inter-filament spaces to likely negligible values (Fuller and Rand 1999; Kabir et al. 2003; Matsuo et al. 2016; Oshima et al. 2016) (Fig. 4A). Depending on values of cross-linker lengths and the presence of hydration shells (or not), bundle mechanical properties may change in consequence, as bundle resistance to stress and buckling, as shown in Mogilner and Rubinstein (2005).

Potentially, free inter-filament spaces will be filled with fluid contained within the filopodium. We made the hypothesis that the cytosol trapped in these inter-filament spaces is taken back to the cell body by the whole bundle. This would further contribute to the outflow volume, thus triggering a correspondingly larger compensatory inflow of matter (cytoplasm) to satisfy the mass conservation law. To test this hypothesis, we need to calculate the flow rate of cytoplasmic fluid entering and leaving filopodium as a function of the inter-filament spaces in the possible configurations that may occur (i.e., with and without hydration). By configurations, we mean bundles with a different number of filaments and length of cross-linkers. Because of the hexagonal packing of filaments in the bundle, a filopodial cross-section area is essential to identify (and quantify) spaces where the cytosol fluid can be stored when outcoming from filopodium.

Let us name $S_{{{\text{empty}}}}$ the fraction of filopodial cross-section area corresponding to inter-filament spaces. Obviously, that depends on the number of filaments N and crosslinks with length s. Being $S_{C}$ the total area corresponding to the circles in Fig. 5 and $S_{H}$ the total area of packed hexagons, we can obtain the packing density $\varphi =\varphi \left(r,s\right)={S}_{C}/{S}_{H}$ covered by the filaments. To do this, we treat two cases for the length of the links: $s = 0$ in case of hydrated filaments, and $s > 0$ in intermediate scenarios or absence of hydration.

Fig. 5

(Available at www.geogebra.org)

Hexagonal packing of actin filaments. In a, a cross-section representation of packed hexagons surrounding filaments with presence of hydration. In b, a cross-section representation of a filopodium in absence of hydration, with the dashed circle illustrating the plasma membrane and a centered linked filament bundle, as a geometric model for the hexagonal pattern exhibited in Fig. 2. In both, the assumed hexagons in the packing are represented to make feasible the calculation of the total cross-sectional area considered in the model as inter-filament spaces, which is highlighted in yellow. The fraction of the packed hexagons’ areas discarded is considered filled by inflow cytoplasm since it is sheltered between the most external hexagons of the packing and the membrane. Proportions are out of scale. Image created by Thiago F. Leal with Geogebra

For $s = 0$, density $\varphi$ does not depend on filament radius r or any other parameter, providing 90.69% occupancy of the packing space. It means that we have 9.31% of “empty area”, where hypothetically the fluid could flow through it. For $s > 0$, that is, in intermediate hydration stage or in its absence, the parameters r and s must be carefully observed because their variation will directly impact all other data.

From there, we have the percentage of regions in the hexagonal packing occupied by filaments in each case. The area of empty regions is easily obtained by doing $\left(1-\varphi \right)$. Then, the total cross-sectional area that represents the inter-filament spaces in a bundle of actin filaments is:

$$ S_{{{\text{empty}}}} = \left( {1 - \varphi } \right) \cdot S_{H} \cdot N - \varepsilon . $$

(4)

Here, we are using ε to represent the fraction of this calculation regarding “free inter-filament spaces” in the bundle. These regions have their areas calculated but do not correspond to what we want. We consider that cytosol can outflow of filopodia just through the inner of the actin filaments bundle. To make it clear, let define a “bundle edge” to delimitate the region that we can consider within the bundle. That boundary joins the centers of most external assumed hexagons in the pack and it is represented by a blue line in Fig. 5B. With this said, to estimate ε we observed some optimal cases of hexagonal packing of circles and the areas out of the packing edge in each case, to better fit in the biological context.

A weighted average between a set of optimal configurations can give us a good hint to obtain ε. Since N is the number of packed hexagons, we can use $N=\{\mathrm{10,19,29,51}\}$, which provides us the respective set of percentual values ${\varepsilon }_{N}=\{\mathrm{40,30,27,20}\}$. We see that, as N increases, ε decreases. But in the biological scenario in question, the number of filaments in a bundle will not much exceed the greater N considered above. In effect, in vivo, they rarely reach such values. Then, it follows:

$$ \varepsilon = \frac{{\sum {N \cdot \varepsilon_{N} } }}{\sum N } = 0.2544 \sim 0.25. $$

From this point onward, we will estimate that about 25% is a disposable area in the calculation of $S_{{{\text{empty}}}}$. In practice, we rewrite $S_{{{\text{empty}}}} =$ $0.75\cdot \left(1-\varphi \right)\cdot {S}_{H}\cdot N$. With these estimates taking place, we begin to determine areas and flow rates that correspond to incoming and outcoming of matter in the filopodium. This is fundamental to quantify the volume flux within filopodium, which causes a drop of pressure due to the volume removal in filopodial tip, thus inducing a compensatory inflow of cytoplasm as a response.

Naming $S_{{{\text{filop}}}}$ the whole cross-sectional area of the filopodium and $S_{{{\text{bundle}}}}$ the total area of the N-filament hexagonal packing, we can calculate the area $S_{in}$ through which cytosol can flow in direction of the tip and write $S_{{{\text{in}}}} = S_{{{\text{filop}}}} - S_{{{\text{bundle}}}}$. From expression (4), we have $S_{H} N - \varphi \;S_{H} N - \varepsilon = S_{{{\text{empty}}}}$. Noting that $\varphi = {{S_{C} } \mathord{\left/ {\vphantom {{S_{C} } {S_{H} }}} \right. \kern-0pt} {S_{H} }}$, it follows that $\left({S}_{H}-{S}_{C}\right)\cdot N-\varepsilon ={S}_{{\text{empty}}}$.

If we compute $\left({S}_{H}-{S}_{C}\right)\cdot N-\varepsilon$, we will get the empty spaces inside the “bundle edge”, as we agreed previously. That will be exactly the way through which the fluid would outflow of filopodium, according to our hypothesis. So, calling that as $S_{{{\text{out}}}}$, we finally have $S_{{{\text{out}}}} = S_{{{\text{empty}}}}$, which we are able to calculate already. Due to the principle of mass conservation, we have continuity equation for fluids, which gives that in a confined system and at a time interval Δt, $Q_{{{\text{in}}}} = Q_{{{\text{out}}}}$ holds, where $Q_{{{\text{in}}}}$ and $Q_{{{\text{out}}}}$ are filopodial inflow and outflow rates, respectively.

The flow velocity is decomposed into parts related to the contribution of each phenomenon. And it is worth noting that each kind of transport mechanism will affect each part of the mixture in a specific way. Thus, on the incoming flow velocity $v_{{{\text{in}}}}$ there is contribution of diffusion and advection processes, with $v_{{{\text{dif}}}}$ velocity due diffusion and $v_{{{\text{drift}}}}$ velocity due fluid dragging.

So, inflow rate may be described as:

$$ Q_{{{\text{in}}}} = S_{{{\text{in}}}} \cdot v_{{{\text{in}}}} ,\;{\text{where}}\;v_{{{\text{in}}}} = \;\;v_{{{\text{dif}}}} + v_{{{\text{drift}}}} . $$

(5)

Similarly, outflow rate will be:

$$ Q_{{{\text{out}}}} = S_{{{\text{out}}}} \cdot v_{{{\text{out}}}} ,\;{\text{where}}\;v_{{{\text{out}}}} = \;\;v_{{{\text{mol}}}} + v_{{{\text{cytosol}}}} . $$

(6)

In this case, outcoming flow velocity $v_{{{\text{out}}}}$ combines velocities of molecules $v_{{{\text{mol}}}}$ (actins and ABPs) and $v_{{{\text{cytosol}}}}$ velocity on which the cytosol, by our hypothesis, outflows through inter-filament spaces, or gets trapped into hydration coats, or both. Molecules will be subjected to actin dynamics and polymerization. Since polymerization maintains F-actin retrograde flow $v_{{{\text{ret}}}} = 70$ nm/s, we make $v_{{{\text{mol}}}} = v_{{{\text{ret}}}}$ for G-actins.

This flux is induced by the generation of a negative pressure in filopodial tip as long as the fluid reorganization occurs through polymerization, where molecules and cytosol are consumed through specific dynamics. The total outflow volume (molecules and trapped cytosol) is precisely computed from the assumptions that it had just been done to provide an overview of the cytoplasmic flux within the filopodium through volumes and velocities in each model phase. Further, the model description must take into account molecular interactions and how filopodia grow as decisive information about how diffusion and drift can be associated in guiding filopodial dynamics.

2.4 The motion engine

2.4.1 Coupling advection and diffusion

Thus, we can start making assumptions about molecules’ movement. Particles of a specie p diffuse within the growing filopodium on a continuous random walk with diffusion coefficient D_p. Despite the chaotic behavior of the motion process, there is a preference for walking along the filopodia axis, but also a mechanism to make possible low and random “up and down” fluctuations in the particle’s trajectory since collisions are not explicitly simulated.

Here, as already said, particles will always move, but react just when they get close enough. It can be seen as an adaptation of the model described in Erban et al. (2014), where each particle moves or reacts, not both, at each time step. The reactions will occur respecting parameters of chemical affinity held in every entity of the system, being particles or filaments. So, consider a particle of type p diffusing in one dimension with diffusion coefficient D_p. According to Islam (2004), the mean time that p takes to diffuse from $x = 0$ to $x = h$ is

$$ \left. {T_{h} } \right|_{p} = \frac{{h^{2} }}{{2D_{p} }}. $$

(7)

If we compute simply $n \cdot \left. {T_{h} } \right|_{p}$, with $n \in {\mathbb{Z}}_{ + }$, then we have the time that the same particle diffuses at a distance $\ell$ = n·h, which brings

$$ \left. {T_{{\ell}} } \right|_{p} = n \cdot \frac{{h^{2} }}{{2D_{p} }} = \frac{{{\ell} \cdot h}}{{2D_{p} }}. $$

(8)

From that, the average diffusion velocity $\left. {v_{{{\text{dif}}}} } \right|_{p}$ of a particle p is

$$ \left. {v_{{{\text{dif}}}} } \right|_{p} = \frac{{\ell}}{{\left. {T_{{\ell}} } \right|_{p} }} = \frac{nh}{{\frac{{nh^{2} }}{{2D_{p} }}}} = \frac{{2D_{p} }}{h}, $$

(9)

which will be treated as a diffusive contribution in the velocity of every particle p. This is in agreement with (Islam 2004) which estimates the approximate time required for a particle to diffuse over a given distance x, in an environment where its diffusion coefficient is D, as $t \approx x^{2} /q_{i} D$, where q_i is 2, 4, or 6 depending on the number of dimensions (i = 1, 2 or 3). In this model, we consider a linear displacement (which states that q₁ = 2) and a G-actin monomer diffusion coefficient D_p = 5 µm²/s. In further simulations with Actin-Profilin dimers, we consider D_p = 6.67 µm²/s in response to their molecular mass.

Coupled to $\left. {v_{{{\text{dif}}}} } \right|_{p}$, we have the contribution of advection to the speed of particles. It is worth remembering that $\left. {v_{{{\text{drift}}}} } \right|_{p}$ is closely related to the flow of volume out of the filopodium, quantified from inter-filament spacing estimate. According to Mogilner and Rubinstein (2005), the drift part of G-actin flux is equal to the speed of filopodium extension, that is:

$$ \frac{{{\text{d}}L}}{{{\text{d}}t}} = \left. {v_{{{\text{drift}}}} } \right|_{{{\text{actin}}}} , $$

(10)

which in other terms can be written as:

$$ \left. {v_{{{\text{drift}}}} } \right|_{{{\text{actin}}}} = k_{{{\text{on}}}} \cdot \delta \cdot C(L) \cdot e^{{\left( {{\raise0.7ex\hbox{${ - N_{0} }$} \!\mathord{\left/ {\vphantom {{ - N_{0} } N}}\right.\kern-0pt} \!\lower0.7ex\hbox{$N$}}} \right)}} . $$

(11)

This is due to the fact that membrane displacement during elongation implies an increase of internal volume that must be filled with input (and drag) of matter. Then, the right term of expression (11) quantifies how much filopodium dynamically increases in length, based on the number of polymerized subunits. Each molecule contributes with $\delta$ in filament length. That depends on the assembly rate $k_{{{\text{on}}}} = 1{0}{\text{.8}}\;\mu {\text{M}}^{ - 1} {\text{s}}^{{ - 1}}$, which in its turn is proportional to the local G-actin concentration.

The exponential factor is responsible for slowing protrusion rate with respect to the number N of filaments, which would require more or less polymerized monomers to grow. Important to notice that $N_{0} /N$ is also related to membrane resistance as a limiting factor of growth, since if $N_{0} > N$, bundle is not enough stiff to perform mechanical effect on membrane. However, when filopodium is in a steady state (with no elongation or retraction), there is still the consumption of matter that induces an incoming flux of cytoplasm by hydrostatic pressure. For such situations, we use information about volume and speed of the outflowing matter and the inflowing flux of cytoplasm for mass conservation purposes.

Transcribing that in the notation used here, it follows:

$$ \left. {v_{{{\text{drift}}}} } \right|_{p} = \left( {\alpha \cdot k_{{{\text{on}}}} \cdot \delta \cdot C(L) \cdot e^{{\left( {{\raise0.7ex\hbox{${ - N_{0} }$} \!\mathord{\left/ {\vphantom {{ - N_{0} } N}}\right.\kern-0pt} \!\lower0.7ex\hbox{$N$}}} \right)}} + \beta \cdot \frac{{S_{{{\text{out}}}} \cdot v_{{{\text{out}}}} }}{{S_{{{\text{in}}}} }}} \right) \cdot M_{p}^{ - 1} . $$

(12)

2.4.2 Computing particles-velocities

In expression (12), drift by elongation and/or advection by flux are considered as possible contributions to particles’ velocities. The parameters α and β are binary and used to control this process. Once there is elongation, $\alpha = 1$; otherwise, $\alpha = 0$. Also, basically $\beta = 1$, but we have the possibility to change it for purposes of virtual experimentation. Then, no physical effect on fluid motion will be lost. Moreover, representing molecular mass of particle p, a factor M_p was set to make speed rates proportional to particle mass. Since profilin, Ena/VASP, and actin do not have same size, it can induce different flux effects on each type of particle.

From that, we define:

$$ \left. {v_{x} } \right|_{p} = \left. {v_{{{\text{dif}}}} } \right|_{p} + \left. {v_{{{\text{drift}}}} } \right|_{p} , $$

(13)

in x-axis direction for every particle p.

Besides, particles may experience possible and random “up and down” fluctuations with a velocity $\left. {v_{y} } \right|_{p}$, and then the velocity vector of particle p may be written as

$$ \overrightarrow {{\left. v \right|_{p} }} = \left. {v_{x} } \right|_{p} \cdot \overrightarrow {{\varvec{i}}} + \left. {v_{y} } \right|_{p} \cdot \overrightarrow {{\varvec{j}}} . $$

(14)

Every mentioned expression here was presented in their dimensionless form. Variables and parameters were properly rescaled due to compartment length and time scale involved during simulations.

So far, it has been established: means for the movement of particles, bases to justify the coupling of the described transport phenomena, foundation to calculate geometric parameters of F-actin bundle, and the systemic description of the reactions involved. Then, the construction of a mathematical formalization through a stochastic bias to describe the problem was presented and it will have its computational approach elaborated in the next section.

3 Fluid flow in filopodia: re-organization of cytosol

In the biological context, the structures involved hold fundamental properties that impact flow dynamics and change alongside. As an example, G-actins react with actin binding proteins (ABP) (Pak et al. 2008; Prokop et al. 2012) following specific chemical affinities and each bond formed among them evolves according to parameters such as reaction’s strength and duration. Furthermore, molecules’ position change stochastically and are influenced by (and may cause) chemical interactions. Likewise, each filament is constituted by a number of actin monomers, whose variation change dynamically mechanical properties of filopodia. Much of filopodia’ architecture and dynamics depend on aspects of chemical bonds whose properties vary with the phenomenon’s behavior. We address these challenges using object oriented programming (OOP) theory and methodology as a modeling paradigm to simulate both molecules, filaments, and chemical bindings.

The representative scheme in Fig. 6 illustrates, in a simplified way, the main interactions between the chemical objects modeled. The stochastic aspect of simulations is present in the particles’ motion towards the filopodial tip and fluid reorganization as posed in Sect. 2. In Fig. 6, circles were used to designate particles (color description follows in the figure caption), which are mostly molecules and squares represent the reactions between any two molecules.

Fig. 6

OOP scheme of actin dynamics in filopodia. A G-actins (orange) have attraction to F-actin tip; B Ena/VASP (blue) binds F-actin tip to prevent capping, while G-actins are sequestered by profilin (green); C as hetero-dimer, profilin has high affinity for Ena/VASP at the tip, promoting G-actin polymerization through a catalytic process; D Once G-actin is polymerized, Ena/VASP jumps to the new F-actin at the tip and profilin is released; E Ena/VASP also has a role in bundling filament tips; F Fascin (yellow) as cross-links actin filaments away from the tip, and ATP-Actin (circle with grey triangle) becomes ADP-Actin (circle with pink triangle); the bundle produces force against the membrane so that filament treadmilling occurs. Images by A. Prokop

In the model proposed here, biological structures are computationally characterized as classes and their objects to simulate all biochemical processes that may impact the filopodial dynamics. Reactions of any kind considered here invariably aim to carry out transport of G-actin and/or favoring polymerization. We will describe only those involving Profilin, Ena/VASP, and actin in globular or filamentous states, but it should be highlighted that more players may act in filopodial dynamics and be modelled in this framework.

We define:

• i. molecule class, with three objects: Actin, Profilin, and Ena/VASP; each particle memorizes their own position (updated at each time step), a label to control how many molecules of each type are present in the simulation, and the attribute “status” to inform if a particle is bound to another (“off”) or free to react (“on”);

• ii. molecular-bond class, which basically means the bindings among particles; objects are formed depending on which particles are involved and to their chemical affinities; possible interactions are polymerization (AA), sequestration by profilin (AP), catalytic reaction between actin-profilin complexes to any Ena/VASP present at the bundle tip favoring actin polymerization (AE then AA);

• iii. bundle class, whose objects are the F-actins resulting of molecule-objects reorganization due to polymerization; F-actins are structured by bipartite graphs for containing molecule and binding objects; polymerization and depolymerization are computationally simulated by inserting or removing data in filament-objects corresponding graphs.

The re-organization of cytosol is a consequence of the polymerization that occurs at the tip of a filopodium (Fig. 3, Sect. 2.2). Our model implements the polymerization process and the consequent organizational phase-transition (see Fig. 3, Sect. 2.2) by considering each molecule suspended in the cytosol as a computational object possessing affinity and other chemical attributes which may vary, besides its location and mass. Object orientation, specifically its run-time features, is a key modelling artefact here. Each inflow molecule suspended in the cytosol, as well as any possible binding, is a computational object in the programming sense. They attach to each other during run-time by the standard method activation of object-oriented programming languages.

Computational molecule objects and chemical bonds are defined as object classes in the program and instantiated whenever something happens, when a molecule enters the domain at its boundary or when a chemical reaction occurs, for instance. Taking into account that molecules have localization attributes, their ‘localized’ run-time configuration represents neatly the fluid organization. Changes in the run-time configuration, induced by chemical reactions and the presence of chemical-bond objects among molecules, represent re-organization. The present version of the model deploys, though, only molecules and biological entities that are strictly necessary to make the flow within filopodia happen. One aspect of the bundle movement not yet represented but which can be easily addressed with this modelling methodology is the dispersion of F-actins at the base of the filopodia into the central cell compartment, where the binding strength plays a crucial role.

Polymerization and depolymerization are computationally simulated by inserting or removing data in filament-objects corresponding graphs as illustrated in the codes below. For instance, the code (4) shows that polymerization of one subunit increases F-actin length in 2.7 nm. On the other hand, depolymerization occurs as a function of the simulation time steps, whose rate is $k_{{{\text{off}}}} = \, 1$ subunit/step on a set of $N_{{{\text{step}}}}$ filaments, what can be seen in code (5).

Code (4):

Code (5):

In the literature, some studies use rates to determine the probability of reactions to occur (Erban et al. 2014). In this present work, reactions will be subjected to proximity among particles, simulating what seems to be reasonable to happen in vivo; two particles with mutual chemical affinities react when they are close enough. The probability of occurring such reactions will depend on the molecules’ own random movement once they flow into the filopodium. Then, position attribute information stored on every particle-objects (defined in the molecules class) are used to determine if two particles must react at some time step.

The code (6) below indicates all procedures described in previous sections and shows the sequence of computational processes in simulations. Initially, classes and objects are instantiated with respective attributes and relevant methods. Once executing the program, the function “Initialize classes” indicates an initial construction of the most basic data structures that will serve to store objects. Two lists are created: one to record all the molecules in the system (inserted by “Create molecules” procedure); another to store information of each filament in the bundle.

Code (6):

We emphasize that the modeling strategy here proposed allows the insertion of new players (such as new ABPs and the consequent reactions that arise) by simply defining new class instances. This favors the analysis of which effects these insertions lead on the entire system, or to test hypotheses about the role of each new introduced object. With this approach, we can simulate filopodial growth dynamics based on a specific set of interacting elements. In future developments, the complexity of individual F-actin simulation can be increased by inserting F-actin cross-linking proteins to observe bundle response to stress and membrane load, as proposed in Peskin et al. (1993), but still keeping the assumption of trapped fluid and hydration shells due to the hexagonal packed pattern, as well as considering the filopodial elongation promoted by the coupling of diffusion and drift as transport mechanisms. This can be seen in the following section.

4 Virtual biological experiments

Below we describe some problems for which the model has proven useful, providing examples of applications that can be redesigned, extended, or improved. The test in Sect. 4.1 exploits the influence of the inter-filament spacing parameter, and variations in the number of bundle filaments to inspect whether they could induce more or less speed in the flow, affecting the elongation of the filopodium. In Sect. 4.2, we tested the predominance of each coupled transport method, as an addendum to what had been predicted in Moura et al. (2016). In Sect. 4.3, we used the model to test hypotheses about the behavior of the output flow combined with obtained data in the two previous tests, which revealed that the cytosol + filament encapsulation thesis is viable. Finally, Sect. 4.4 is dedicated to exhibit one of the reasons for choosing OOP as a modeling paradigm; namely, the possibility of easily making different type of particles in the system interact, inserting new ones, or isolating them to verify their influence and effect on the dynamics.

4.1 Inter-filament spacing and filopodial growth

A sample of values for N = {20, 30, 50} was chosen to estimate the variation of the filopodial length L (in micrometers) when inter-filament spacing is changed from $s = 0$ to $s = 12$(in nanometers). Two fundamental scenarios have been considered: first, the fluid mixture contains only G-actins; second, actin monomers are regulated by ABPs. At this point, active transport by motors is not involved.

Figure 7 shows the distributions of 1500 and 3000 G-actin particles at the end of a default time of simulation ts = 100. The distributed particles in red refer to moving G-actins. The domain represents a filopodial tube with length in micrometers and diameter in nanometers. F-actin bundle has N filaments computed individually, but it is embodied on a black rod reflecting the bundle actual length. To improve visualization, the distance between the bundle tip and membrane leading edge was chosen to be set as 0.5 µm. The axes are not in their natural ratios.

Fig. 7

G-actins distribution combining values of N = {20, 30, 50}, s = {0, 12}, and N_P = {1500, 3000}. Lengths are those in Table 1 and the axes were re-scaled for comparison

Table 1 Variation of L as function of s

Two other points must be highlighted. First, the increase of $N_{P}$ leads to higher values of L for all N. This can be physically explained by a greater supply of G-actins to polymerization. Second, for each fixed number of filaments, the filopodial length L increases whenever inter-filament spacing changes from $s = 0$ to $s = 12$. This bears a physical interpretation: an increase in outflowing volume induces an increase in the number of drifting particles. This occurs since the cross-section outflowing area $S_{{{\text{out}}}}$ increases along changes in s. In both points, higher removal of volume provokes a decrease in the outflow velocity, which produces as a consequence lower rate of filopodial elongation.

To confirm those facts, we must analyze how filopodial dynamics behaves in response to the biological regulation of ABPs. In Fig. 8, the yellow particles refer to moving hetero-dimers AP, in amounts of 1500 and 3000 particles. The formation of hetero-dimers throughout filopodial interior is not taken into account; the inflowing particles are already dimerized.

Fig. 8

Hetero-dimers AP distribution combining values of N = {20, 30, 50}, s = {0, 12}, and N_P = {1500, 3000}. Lengths are those in Table 2 and the axes were re-scaled for comparison

From Tables 1 and 2, interesting knowledge should be gained or highlighted. The main contributions of this section are:

• • the verification of an inverse relation between the number of filaments N and filopodial length L, already approached by Mogilner and Rubinstein (2005) although focusing on protrusions and diffusion as the transport mode;

• • a novel relation between filopodial length L and inter-filament spacing s, where L increases whenever s changes from 0 to 12; although counterintuitive, as a greater space for the cytosol to outflow involves a smaller flux speed, but does not affect the proper supply of polymerization;

• • the suggestion of a new hypothesis to be investigated, which is the direct relation between filopodial length L and the number of moving particles N_P; the increase of N_P enhances filopodial growth in all tested scenario, despite a reduction in the percentage of length variation with respect to s.

Table 2 Variation of L as function of s in hetero-dimers AP distributions

4.1.1 Statistical evaluation of model stability

In the previous sub-section, the presented results were chosen from several performed virtual experiments to illustrate particles distributions. It should be clear that new runs generate different distributions mainly due to the randomness of the particles’ movement. Besides, the pattern of results is not much beyond those shown, with average filopodial lengths not varying more than ± 0.1 µm. To probe the stability of the model, Table 3 shows the variance σ between the final lengths of filopodia after 10 runs for each set of parameters tested in Sect. 4.1. As the methodology of variance calculation is maintained for all cases, the evaluation of scenarios with G-actins and ABPs will be omitted.

Table 3 Variance between the lengths after 10 runs for each considered case

4.2 Predominant transport phenomenon: diffusion vs advection

The spatial distribution of particles within filopodia provides a picture of how important inter-filament spacing are, the number of particles, and time in filopodial growth when diffusion, drift, and those reactions relevant to actin dynamics are considered. Nevertheless, the mechanisms that drive particle flow along the filopodium need to be investigated. An immediate question in this direction is the effect of each phenomenon separately.

In the following, a broader understanding of the transport mechanisms is provided by inspecting velocity profiles of the inflow mixture for given conditions. This allows us to understand which of the three processes is more effective and where. We also analyze their individual contributions and how these transport phenomena relate to fundamental parameters such as the inter-filament spacing s, the number of filaments in the bundle N, the length of the filopodium L, or chemical aspects like biological regulation.

The following curves represent velocity profiles of inflowing cytoplasm with G-actins dispersed. Figures 9, 10, and 11 display together the curves for different N values for each mode: diffusion, drift, and inflow velocity, with $N = \{ 20,\;30,\;50\}$ and s = 0. The x-axis registers the filopodial length, in micrometers. Flow velocities given by each phenomenon along filopodial length are set on y-axis, in micrometers per second. The total simulation time for all profiles is t_s = 100.

Fig. 9

Diffusion contribution to flow velocity in a filopodium with only G-actins. As long as filopodium elongates, diffusion rapidly becomes inefficient to sustain G-actin transport. Also, greater number of filaments diminishes filopodial final length

Fig. 10

Drift contribution to flow velocity in a filopodium with only G-actins. As long as filopodium elongates, with consequent diffusion decay, drift becomes the main physical transport phenomenon, with higher velocity for lower number of filaments

Fig. 11

Total velocity of inflow cytoplasm in a filopodium with only G-actins. As long as filopodium elongates, with consequent diffusion decay, drift becomes the prominent physical transport phenomenon, with higher velocity for lower number of filaments

These figures point to important features of the physical phenomena involving particle movements in filopodia. First (Fig. 9), diffusion alone quickly becomes fairly ineffective to deliver actin to polymerization. The speed absolute value tends to zero once filopodial length becomes higher than 1 µm for all three values of N. Further, the relation between N and L discussed in Sect. 4.1 is perceived: length becomes longer for lower number of filaments in the bundle, and vice versa.

The fast decay of diffusion velocity is a known result (Mogilner and Rubinstein 2005), but it is interesting to see that diffusion is related to filopodial configuration. An important question is how diffusion relates to drift, which has its contribution for chosen N deployed in Fig. 10. In the virtual experiments, drift increases for almost the entire process, except for a short time in the initial stage of growth when its value is not significant in comparison to diffusion.

In Fig. 11, we saw the composition of both phenomena, i.e., diffusion and drift. This composition means the total inflow cytoplasm velocity $v_{{{\text{in}}}}$, defined in Sect. 2.2 as $v_{{{\text{in}}}} = v_{{{\text{dif}}}} + v_{{{\text{drift}}}}$. Besides, the total speed $v_{in}$ can be interpreted by two central features of the filopodial domain: the source of matter at the filopodial base and the sink in the polymerization point, which induces an increasing advective effect due to the volume removal and consequent drop of hydrodynamic pressure at filopodial tip. Hence, the curve $v_{{{\text{in}}}}$ is approximated by diffusion curve in the beginning of the filopodial tube, while $v_{{{\text{in}}}}$ approaches drift magnitude in the end of it. For each value of N, drift becomes more prominent as diffusion value tends to zero.

These patterns emerge for all parameter sets, with greater or smaller contribution of different mechanisms at specific locations. For example, for N = 30, the drop of diffusion curve and the rise of drift curve are slightly more pronounced, which is caused by factors as: a higher demand of G-actins by polymerization, a greater volume of hydrated cytosol, and/or a larger inter-filament spacing for the cytosol to outflow. In comparison, due to the inverse relation between N and L, the total speed leads to similar dynamics for N = 20, but with its curve exhibiting lower growth close to the sink.

The velocity profiles for s = 12 are similar to those with s = 0, but a perceptible difference is the slightly higher intensity of drift for all N values, especially N = 50. This is related to the discussion in the previous section about greater supply for polymerization providing greater consumption of matter, which induces more drag of cytoplasm. At last, analogous observations can be done when ABPs are present in the domain and curves display patterns deeply similar to those in the scenarios analyzed individually, with solely G-actins. The most relevant contribution of ABPs lays in terms of larger maximum filopodial lengths when considering the default time of observation t_s = 100.

4.3 Encapsulated cytosol outflow is a feasible explanation

In the present investigation, we establish a hypothesis for cases of total hydration, that is, $s = 0$. Imagine that the entire bundle is filled with water (colloid) associated with the filaments by hydration. In this situation, there would be no flow through the bundle spaces. Every volume that flows out of the filopodium is given by the filaments, and cytosol as a hydration shell that occupies the entire interior of the bundle and F-actin surroundings.

For this, we consider that the encapsulated outflowing volume would be given by

$$ V_{{{\text{outflow}}}} = V_{{{\text{bundle}}}} = 0.75 \cdot S_{H} \cdot N \cdot L_{c} , $$

(15)

and outflows at the speed of the retrograde flow $v_{ret}$. The calculations aim to answer whether this dynamic is feasible in these terms.

The initial conditions used in this test are: (A) only G-actins in the system, (B) presence of ABPs, and (C) G-actins, ABPs, and active transport by motors.¹ For each situation, the computations seek to answer three questions:

(1) Is incoming volume of cytosol (mixture) equal to volume of outflowing encapsulated bundle with hydration shells?

Answering this question shows if inflow and outflow volumes are in balance, satisfying mass conservation. The answer is “yes” just if the volumes are balanced during all the simulation time. A tolerance of 0.1 was set for the difference between those values to prevent approximation errors.

(2) Is polymerization properly supplied by actin monomers?

This simulation calculates the percentage of cases where $N_{{\text{p}}} > N_{{{\text{ret}}}}$ (where $N_{{{\text{ret}}}}$ is the number of actins needed to feed F-actin retrograde flux, as in Mogilner and Rubinstein (2005)). The answer is “yes” if polymerization is supplied during all the simulation time.

(3) What is inflowing flux velocity?

Parameter to be compared with the results of other cases.

For simplicity, in this first approximation, we do not consider water molecules “bind-unbind” reactions with F-actins at any point of its length. As polymerization, hydration occurs at filaments’ tips. But since the initial filopodial length is zero, the bundle will be hydrated during the whole process. It was chosen N = 30 to perform this simulation.

Then, the results of the computation were the following:

	Question 1	Question 2	Question 3
Scenario A	Yes (within the tolerance)	No % of “yes”: 60.8 % of “no”: 39.2	$v_{{{\text{in}}}} = 0.509\;\upmu {\text{m/s}}$
Scenario B	Yes (within the tolerance)	No % of “yes”: 81.5 % of “no”: 17.5	$v_{{{\text{in}}}} = 0.677\;\upmu {\text{m/s}}$
Scenario C	Yes (within the tolerance)	Yes % of “yes”: 100.0	$v_{{{\text{in}}}} = 0.83\;\upmu {\text{m/s}}$

The question 2 of the last test brings biologically relevant information that leads to further investigation. For the case where all biological interesting dynamic aspects are considered, we verify that the hydration of filaments may not be simply an osmotic cell control mechanism, but also sustains the flow of matter in filopodia. Hydration shells around filaments displacing together with the F-actins treadmilling process contribute to keep the balance of mass within the filopodium.

This assumption does not contradict the previous hypotheses and results, as well as the vast field of knowledge available in the literature. It gives us a perspective of extreme cytosol flow patterns—with or without total hydration. When the spaces in the bundle are larger, the total volume of outflow cytosol can be decomposed into a colloidal flux (whose speed tends to be much greater than $v_{{{\text{ret}}}}$) and hydration shells around filaments (flowing with speed $v_{{{\text{ret}}}}$). This further increases the inflow of molecules to support polymerization.

This is a novel and relevant finding, which paves the way for further mind-boggling investigations. This picture can be widely improved by simulating a completely particulate fluid representing colloidal water molecules flowing with charged actins and their regulators. In this case, the formation of the hydration shells can be simulated and to enhance the understanding of the previous analyzes proposed.

4.4 Influence of ABPs in filopodial growth

The purpose of this application is to evaluate specific scenarios of filopodial dynamics, where particular reacting particles are “turned on” or “turned off” to observe differences in system behavior. This leads to a deeper understanding of the influence of the presence or absence of several combinations of actin-regulator molecules and by observing how the filopodium evolves in each case. It is important to notice that some of these results have been discussed in the literature (Bear and Gertler 2009; Birbach 2008; Harker et al. 2019), but these simulations also demonstrate that our model is a viable tool to investigate these and similar biological interdependencies. For reference, the red line in the following graphs represents $L_{{{\text{crit}}}}$, which is the length where the monomer concentration becomes critical and polymerization rate is not supplied efficiently. Then, in the following tests, we investigate profilin and Ena/VASP influences on filopodial growth by including ("on") or eliminating ("off") them from the system. We choose N = 30, and s = 12 for all tests.

• Basic filopodial model (Ena/VASP = ‘off’ and Profilin = ‘off’): the main ABPs are not present, and G-actins are responsible for the dynamics influenced solely by transport physical processes. See Fig. 12.

• Filopodial growth is improved (Ena/VASP = ‘on’ and Profilin = ‘on’): both regulators taking place; filopodial elongation is improved. See Fig. 13.

• G-actin is sequestered (Ena/VASP = ‘off’ and Profilin = ‘on’): G-actin sequestration decreases polymerization rate, which shortens filopodial length. See Fig. 14 for different tested amount of profilin.

• Decrease in polymerization rate (Ena/VASP = ‘on’ and Profilin = ‘off’): as Ena/VASP binds F-actin tips, the polymerization rate decreases drastically, since no profilin are present to interact with Ena/VASP. But when there is a number of free G-actins, polymerization can still occur. See Fig. 15 to observe scenarios with different amounts of profilin.

Fig. 12

Filopodial length with no ABPs in the domain. Simulation with t_s = [0, 200]

Fig. 13

Filopodial length with mediation of profilin and Ena/VASP. Filopodial lengths are highly impacted by the regulating proteins of G-actins. Simulation with t_s = [0, 200]

Fig. 14

Influence of profilin. Four values for percent of actin sequestration due to the AP reactions were tested: in blue, 0% (all the G-actins are free); in green, 20%; in yellow, 50%; and in magenta, 60%. The black line marks $x = 0$. As more G-actins are in association with profilin, less actins are polymerized, which shortens filopodial length. Simulation with t_s = [0,1000]

Fig. 15

Ena/VASP regulating bundle tip with no profilin in the system. In this test, just G-actins are dispersed in the fluid and it is simulated a binding of Ena/VASP in bundle tip at t_S = 100. It causes a rapid reduction of polymerization if we consider that G-actin must be in a hetero-dimer with profilin to react with Ena/VASP and, as consequence, polymerize actin monomers

One major reference about stochastic modeling of filopodial dynamics is (Lan and Papoian 2008) with which the results showed here can be validated and compared. The main difference concerns the strategies to simulate polymerization on each F-actin filament. As we have chosen to equalize them by distributing the monomers preferentially into filaments of shorter length, so that the bundle reaches the same length and, from then on, elongate, in Lan and Papoian (2008), polymerization in each individual filament is affected by the stress that the membrane causes against its growth and the observations are derived from this hypothesis. However, even with unlike approaches, the variation in filopodium length over time is comparable with the results of Sect. 4.4 in respect to the phenomenon itself, since the proposal of this work considered effect of drift coupled to diffusion in G-actins transport, which results in a greater growth of simulated filopodia.

5 Conclusion and future work

Filopodia, like any other living phenomenon, is an open, ever moving, and ever-changing system. They bend and grow in several directions, regulated by processes yet barely known. Inside filopodia, there is nothing but motion and change delimited solely by cell membrane, which changes its geometry moving constantly around. Biologists accept that what happens inside filopodia creates and sustains its form and also supports its deformations and redirections with respect to its surroundings. To foster knowledge about what drives protrusions to one or another direction or what slows-down or speeds-up movements, processes, and changes, a throughout understanding of how processes within filopodia yield to (external) stimuli is necessary. The internal processes are physical and chemical, centered on movement and recombination, that is, re-organization, of molecules. Filopodia mechanics, elongation, retraction, and bending, depend on these internal processes. Gathering knowledge about them and encapsulating this knowledge at a higher aggregation level is of prime importance to focus attention in what really matters and to see filopodia as a responsive and adapting entity, connecting and interacting with surrounding biological entities. Understanding filopodia and the cytoskeleton at large is an important biological and medical quest. Modelling is a cornerstone to achieve this.

Internally, filopodia organization and processes depend crucially in actin molecules and their polymerization, which occurs at the tip of the filopodia. Actins exist in two major conformational states: G-actins and F-actins. The majority of studies about filopodia take for granted that diffusion in the cytoplasm would transport G-actins to its tip, no matter its length, where polymerization transform them into F-actins. Yet, we have shown that diffusion cannot do that for filopodia protrusions greater than 0.5 µm when real ones may reach several times (110×) this length (Moura et al. 2016). Hence, the diffusion phenomenon alone is not sufficient to explain the inner workings of filopodia. In the same article, we proposed a schematic model that could explain greater elongations, grounded on a radically new perspective: re-organization of an ‘almost’ laminar in-flow of cytosol with suspended G-actins by polymerization into an out-going cytosol flow trapped in a F-actin tube that brings matter back to the cell’s main compartment, where the F-actin bundle is disassembled recovering the G-actin-cytosol suspension (Fig. 3, Sect. 2.2).

This text describes a computational version of this hypothetic model, with two main purposes: verifying if the flow would really transport enough G-actins to the tip and to showcase a program that could be effortlessly modified to include the disassembling moment, as well as other filopodia features, serving as a virtual laboratory for filopodial investigations. It addresses the deformation of domain (Sect. 2.1) and three moments of the above cycle, namely, the in-flow (Sect. 2.2), the out-flow (Sect. 2.3) and the re-organization by polymerization (Sect. 3). Movement in the whole cycle, including the disassembling phase in the cell’s main compartment is subject to diffusion–advection kinematics (Sect. 2.4), emphasizing the organization/re-organization nature of this phenomenon.

This “modelling exercise” successfully implements our 2016 hypothesis about the flow within filopodia protrusions (Moura et al. 2016) and computes elongations that are 20 times longer than current diffusion models in the literature are able to do. Moreover, these elongations were only limited by computer resources available and there is nothing in the model precluding them to grow indefinitely. Therefore, this achievement not only proves that the re-organization hypothesis is feasible and improves current descriptions of filopodia flow, but also emphasizes that organization is a feature of living phenomena that should be seriously considered by the community to understand them, as many biologists propose (Harold 2005; Kritz et al. 2010). Moreover, the use of object-orientation as a modelling rather than just programming paradigm (Sect. 3) paves the way for constructing malleable programs that can be readily re-used to implement even richer experiments to support the investigation of the myriad of questions about filopodia and cytoskeleton that still remain open.

Nevertheless, one last challenge remains to render this computational model fully useful in biological enquiries about cell-protrusions, related either to filopodia or other aspects of the cytoskeleton. This challenge does not relate to the physics, chemistry, or biology of these phenomena. Instead, it is a geometrical challenge regarding the computational description and tracing of a swiftly changing and bending volume. We need to develop a form-tracing algorithm for the full volume of the filopodia that can be seamlessly blended with tessellation and particle-in-cell methods. Not needed to emphasize how difficult this challenge is.