Actin droplet machine

The actin droplet machine is a computer model of a three-dimensional network of actin bundles developed in a droplet of a physiological solution, which implements mappings of sets of binary strings. The actin bundle network is conductive to travelling excitations, i.e. impulses. The machine is interfaced with an arbitrary selected set of k electrodes through which stimuli, binary strings of length k represented by impulses generated on the electrodes, are applied and responses are recorded. The responses are recorded in a form of impulses and then converted to binary strings. The machine’s state is a binary string of length k: if there is an impulse recorded on the ith electrode, there is a ‘1’ in the ith position of the string, and ‘0’ otherwise. We present a design of the machine and analyse its state transition graphs. We envisage that actin droplet machines could form an elementary processor of future massive parallel computers made from biopolymers.


Methods
The overall approach is the following: we simulate the actin bundle network using three-dimensional arrays of finite-state machines, cellular automata. We select several domains of the network and assign them as inputs and outputs. We represent Boolean logic values with spikes of electrical activity, which are schematically represented as a virtual experiment in figure 1. We stimulate the network with all possible configurations of input strings and record spikes on the outputs. Based on the mapping of configurations of input spikes to output spikes, we reconstruct logical functions implemented by the network. In our design of the actin droplet machine, we consider outputs recorded on all electrodes at a given time step as a binary string and then represent the actin droplet machine as a finite-state machine whose states are binary strings of a given length.

Three-dimensional actin network
As a template for our actin droplet machine, we used an actual three-dimensional actin bundle network produced in laboratory experiments with purified proteins (figure 2). The underlying experimental method was shown to reliably produce regularly spaced bundle networks from homogeneous filament solutions inside small isolated droplets in the absence of molecular motor-driven processes or other accessory proteins [15]. These structures effectively form very stable and long-living threedimensional networks, which can be readily imaged with confocal microscopy resulting in stacks of optical two-dimensional slices (figure 2). Dimensions of the network are the following: size along x royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 191135 coordinate is 225 μm (width), along y coordinate is 222 μm (height), along z coordinate is 112 μm (depth), voxel width is 0.22 μm, height 0.22 μm and depth 4 μm.
Original image: A z = (a ijz ) 1≤i,j≤n,1≤z≤m , a ijz ∈ {r ijz , g ijz , b ijz }, where n = 1024, m = 30, r ijz , g ijz , b ijz are RGB values of the element at ijz, 1 ≤ r ijz , g ijz , b ijz ≤ 255 was converted to a conductive matrix C = (c ijz ) 1≤i,j≤n,1≤z≤m as follows: c ijz = 1 if r ijz > 40, g ijz > 19 and b ijz > 19. The conductive matrices are shown in figure 3. The three-dimensional conductive matrix is compressed along the z-axis to reduce consumption of computational resources, scenario of the non-compressed matrix will be considered in future papers.

Automaton model
To model activity of an actin bundle network we represent it as an automaton A ¼ hC, Q, r, h, u, di. C , Z 3 is a set of voxels, or a conductive matrix C defined in §2.1. Each voxel p ∈ C takes states from the set Q ¼ {w , †, }, excited (w), refractory ( †), resting (°) and is complemented by a counter h p to handle the temporal decay of the refractory state. Following discrete time steps, each voxel p updates its state depending on its current state and the states of its neighbourhood u( p) = {q ∈ C : d( p, q) ≤ r}, where d( p, q) is an Euclidean distance between voxels p and q; r ∈ N is a neighbourhood radius. θ ∈ N is an excitation threshold and δ ∈ N is refractory delay. All voxels update their states in parallel and by the same rule:  royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 191135 state † at the next time step t + 1 and at the same moment a counter of refractory state h p is set to the refractory delay δ. The counter is decremented, h tþ1 p ¼ h t p À 1 at each iteration until it becomes 0. When the counter h p becomes zero the voxel p returns to the resting state°. For all results shown in this manuscript, the neighbourhood radius was set to r = 3. Choices of θ and δ are considered in §2.4.

Interfacing with the network
To stimulate the network and to record activity of the network, we assigned several domains of C as electrodes. We calculated a potential p t x at an electrode location c ∈ C as p c ¼ jz : d(c, z) , r e and z t ¼wj, where d(c, z) is an Euclidean distance between sites x and z in three-dimensional space. We have chosen an electrode radius of r e = 4 voxels and conducted two families of experiments with two configurations of electrodes.
In the first family of experiments E 1 , we studied frequencies of two-input-one-output Boolean functions implementable in the network. We used 10 electrodes, their coordinates are listed in table 1 and a configuration is shown in figure 4a. Electrodes E 0 representing input x and E 9 representing input y are Figure 2. Exemplary z-slices of a three-dimensional actin bundle network reconstructed as described in [15]. An example of the network spiking activity as a response to stimulation is shown in figure 1. We stimulated the network with all possible configurations of inputs, recorded the network's electrical dynamics and then extracted logical gates as follows. For each possible combination (i, j, k), 1 ≤ i, j, k ≤ 6, i ≠ j, i ≠ k, j ≠ k, we considered electrodes E i and E j to be inputs, representing Boolean variables x and y, respectively, and electrode E k as output electrode, representing the result of a Boolean function.
To input x =TRUE, we applied a current to electrode E i , to input y =TRUE to electrode E j . Then we recorded the potential at electrode E k . Two-input-one-output logical functions were extracted from the spiking events as follows. Assume each spike represents logical TRUE and that spikes being less than six iterations closer to each other happen at the same moment. Then a representation of gates by spikes and their combination will be as shown in figure 6.
For each combination (ρ, θ), we counted the numbers of gates OR (x + y), AND (xy), XOR (xÈ), NOT-AND (xy), AND-NOT (xy) and SELECT (x, y). We found that overall the total number of gates ν(θ) realized by the network decreases with increase of θ (figure 7a). The function ν(θ) is nonlinear and could be adequately described by a five degree polynomial. The function reaches its maximal value at θ = 7 (figure 7a). OR gates are most commonly realized at θ = 11, AND gates at θ = 6 and xor gates at θ = 5 as well as θ = 7   (figure 7b). A number of AND-NOT gates implemented by the network reaches its highest value at θ = 6 then drops sharply after θ 8 (figure 7c). NOT-AND gates are more common at θ = 5, 7, 9, 11, while SELECT(x) has its peak at θ = 7 and SELECT( y) at θ = 8, 9 (figure 7c). A total number of gates realized in the network with the excitability threshold fixed to θ = 7 decreases with the increase of δ. Oscillations   7d). The three highest values of ν(δ) are achieved at δ = 10, 17 and 20. Let us look now at the dependence of the numbers of OR, AND and XOR gates of the refractory delay δ in figure 7e. The number of OR gates increases with δ increasing from 10 to 15, but then drops substantially at δ = 18 to reach its maximum at δ = 19. Numbers of gates AND and XOR behave similarly to each other. They both have a pronounced peak at δ = 20 (figure 7e).
The gate frequency analysis presented in this section allows us to choose θ = 7 and δ = 20 for an actin droplet machine constructed in the next section.

Actin droplet machine
An actin droplet machine is defined as a tuple M ¼ hA, k, E, S, Fi, where A is an actin network automaton, defined in §2.2, k is a number of electrodes, E is a configuration of electrodes, S = {0, 1} k , F is a state-transition function F : S → S that implements a mapping between sets of all possible configurations of binary strings of length k. In the experiments reported here k = 6.
In our experiments, we have chosen six electrodes, their locations are shown in figure 4b and exact coordinates in table 2. Thus, F : {0, 1} 6 → {0, 1} 6 and the machine M has 64 states. We represent the inputs and the machine states in decimal encoding. Spikes detected in response to every input from {0, 1} 6 are shown in figure 8.
Global transition graphs of M for selected inputs are shown in figure 9. Nodes of the graphs are states of M, edges show transitions between the states. These directed graphs are defined as follows. There is an edge from node a to node b if there is such 1 ≤ t ≤ 1000 that M t ¼ a and M tþ1 ¼ b.
Let us now define a weighted global transition graph G ¼ hQ, E, wi, where Q is a set of nodes (isomorphic to the {0, 1} 6 ), and E is a set of edges, and weighting function w : E → [0, 1] assigning a number of a unit interval to each edge. Let a, b ∈ Q and e(a, b) ∈ E then a normalized weight is calculated as w(e(a, b) with χ takes value '1' when the conditions are true and '0' otherwise. In words, w(e (a, b)) is a number of transitions from a to b observed in the evolution of M for all possible inputs from Q during time interval T normalized by a total number of transition from a to all other nodes. The graph G is visualized in figure 10a.  e(a, c))|e(a, c) ∈ E}. That is for each node we select an outgoing edge with maximum weight. The graph is a tree, see   Figure 11b shows a number of different nodes, generated in evolution of M, stimulated by a given input. There are below 15 different states found in the evolution in responses to inputs 1 to 21 (21 corresponds to binary input string 010101); then a number of different nodes stay around 25. The diagram figure 11c shows how many inputs might lead to a given state/node of M. Some of the states/nodes are seen to be Garden-of-Eden configurations E (nodes without predecessors) and thus could not be generated by stimulating M by sequences from Q − E.
Assume T is a set of temporal moments when the machine responded at least to one input string with a non-zero state. Configurations at each transition t can be considered as outputs representing the function g : {0,1} 6 → {0,1} 6 . As we can see in table 3, transitions at t = 41 and t = 53 correspond to the highest number of different binary strings (e 1 , …, e 6 ). The graph corresponding to g(41) at t = 41 is shown in figure 12 and is not connected. The small component consists of fixed point 40 (string '101000') with two leafs 39 ('100111') and 38 ('100110'). The largest component has a tree structure at e 0 : f 0 (x 0 , . . . ,  Table 3. Fifty-four state transitions of M over all possible inputs: t is a transition step, μ(t) is a number of different states appeared over all possible inputs, P(t) is a set of nodes appeared at t.

Discussion
Early concepts of sub-cellular computing on cytoskeleton networks as microtubule automata [44][45][46] and information processing in actin-tubulin networks [47] did not specify what type of 'computation' or 'information processing' the cytoskeleton networks could execute and how exactly they do this. We implemented several concrete implementations of logical gates and functions on a single actin filament [48] and on an intersection of several actin filaments [49] via collisions between solitons. We also used a reservoir-computing-like approach to discover functions on a single actin unit [50] and  filament [51]. Later, we realized that it might be unrealistic to expect someone to initiate and record travelling localizations (solitons, impulses) on a single actin filament. Therefore, we developed a numerical model of spikes propagating on a network of actin filament bundles and demonstrated that such a network can implement Boolean gates [33].
In the present paper, we reconsidered the whole idea of the information processing on actin networks and designed an actin droplet machine. The machine is a model of a three-dimensional network, based on an experimental network developed in a droplet, which executes mapping F of a space of binary strings of length k on itself. The machine acts as a finite state machine, which behaviour at a low level is governed by localizations travelling along the networks and interacting with each other. By focusing on a single element of a string, i.e. a single location of an electrode, we can reconstruct k functions with k arguments, as we have exemplified at the end of the §3. The exact structure of each k-ary function is determined by F, which, in turn, is determined by the exact architecture of a threedimensional actin network and a configuration of electrodes.
Thus, potential future directions could be in detailed analysis of possible architectures of actin networks developed in laboratory experiments and evaluation on how far an exact configuration of electrodes affects the structure of mapping F and corresponding distribution of functions implementable by the actin droplet machine. The ultimate goal would be to implement actin droplet machines in laboratory experiments and to cascade several machines into a multi-processors computing architecture.
Conventional hardware is static. Actin networks reconfigure dynamically: some filaments disappear by depolymerization, new filaments appear by polymerization. This is not a disadvantage of the actin network computers because: (1) they operate with the speed several orders more than actin treadmilling rate, (2) actin networks can be stabilized, (3) we can employ dynamic reconfigurablity in the computation.
A computation in actin bundle networks is implemented with travelling mechanical or electrical signals. Thus, we could estimate the speed of the signals propagation would be 10 6 μm s −1 , for sound solitons, or 10 5 -10 8 μm s −1 , for action potential speed [52]. Let us take the lowest estimate 10 5 μm s −1 . Assuming maximum linear size of an actin droplet machine is ca 250 μm, the machine can process about 400 parallel inputs per second, thus operating at 0.4 kHz frequency. Commonly, actin polymerization speed is estimated to be 4 × 10 −1 μm s −1 [53]. An acting bundle has up to 500 actin filaments, which will not fail simultaneously. In fact, we have seen that the networks, once formed, could remain stable over hours without major rearrangements. In contrast with cells, no actin accessory proteins were used and no energy in the form of ATP was provided. The structures selfassembled solely driven by thermodynamic arguments into a stable, frozen state. If we would neglect these experimental findings and would assume an active network with high treadmilling rates, we can consider the network being fixed for at least 10 s, which allows us to execute up to 4 × 10 3 cycles of computation.
The life-time of the fixed network can be even substantially changed by using accessory proteins such as purely synthetic actin crosslinkers from DNA and peptides [54], increasing a ratio of integrine [55] and drebrin [56] peptides in the matrix solution, hardening the filaments with α-actinin [57] and stabilizing the filament with synthetic mini-nebuline [58,59]. Using accessory proteins such as gelsolin, cofilin, formin and myosins would even allow us to speed up potential reconfiguration effects, enabling us to build up a dynamic computing system [6].
Dynamical reconfiguration of actin network computers can be used as an advantage for accelerating Boolean satisfiability solvers [60], reconfigurable data flow machine for implementing atomic functional programming languages [61], dynamical genetic programming on evolvable Boolean networks [62,63], and cryptographic applications [64].
Data accessibility. Supporting data are available at [43]. Authors' contributions. F.H. and J.S. provided experimental laboratory data. A.A. designed the model. All authors analysed the results and co-wrote the paper.