Hiding the squid: patterns in artificial cephalopod skin

Cephalopods employ their chromomorphic skins for rapid and versatile active camouflage and signalling effects. This is achieved using dense networks of pigmented, muscle-driven chromatophore cells which are neurally stimulated to actuate and affect local skin colouring. This allows cephalopods to adopt numerous dynamic and complex skin patterns, most commonly used to blend into the environment or to communicate with other animals. Our ultimate goal is to create an artificial skin that can mimic such pattern generation techniques, and that could produce a host of novel and compliant devices such as cloaking suits and dynamic illuminated clothing. This paper presents the design, mathematical modelling and analysis of a dynamic biomimetic pattern generation system using bioinspired artificial chromatophores. The artificial skin is made from electroactive dielectric elastomer: a soft, planar-actuating smart material that we show can be effective at mimicking the actuation of biological chromatophores. The proposed system achieves dynamic pattern generation by imposing simple local rules into the artificial chromatophore cells so that they can sense their surroundings in order to manipulate their actuation. By modelling sets of artificial chromatophores in linear arrays of cells, we explore the capability of the system to generate a variety of dynamic pattern types. We show that it is possible to mimic patterning seen in cephalopods, such as the passing cloud display, and other complex dynamic patterning.


Introduction
In this paper, we present an application of smart materials, inspired by biological chromatophores, to generate active dynamic patterns. Chromatophores are small pigment-containing cells embedded into the skin of animals, such as amphibians, reptiles and fish. In particular, the cephalopod employs its skin for controlled chromomorphic effects, serving as effective tools for signalling and camouflage [1]. For example, the Sepia apama is known for its 'passing cloud' display, where bands of blue-green colour propagate as waves across the skin, as shown in figure 1a. This visual effect acts to distract and divert predators. Individual chromatophores actuate under neuro-electrical stimulus of flat wedge-shaped muscles [4], causing the pigmented region (saccule) to increase in area and affect local skin colouring, as illustrated in figure 1b. The tight, neuro-muscular connections coordinate actuation, allowing cephalopods to rapidly shift between skin colours and pattern types.
Biomimetic cephalopod chromomorphism allows for new approaches to camouflage and signalling using soft and compliant artificial skin. Compared with other proposed active camouflage solutions, such as retroreflective projection technology [5], an artificial skin requires no external projectors, making it easier to maintain patterns when mobile. By employing complex and dynamic patterns, users may stand out in times of danger, useful for signalling applications such as search and rescue operations.
One promising class of materials for creating artificial muscle is electroactive polymers (EAPs), smart polymers that exhibit a change in shape or size in response to electrical stimulus. Their properties are well suited to mimicking muscular systems in animals since they can be engineered into lightweight and compliant actuators with relatively high work output [6]. Many  using ionic polymer-metal composites (IPMCs) [7] and the swimming motion of eels using polymer gels [8].
Dielectric elastomers (DEs), a type of planar-actuating EAP, offer potential towards realizing effective artificial cephalopod skin. Previous studies investigating a three-spot artificial DE chromatophore system (shown in figure 1c) indicate significant optical modulation that is of the same order as cephalopod chromatophores [3]. Furthermore, their compliance and rapid actuation speed should allow for wearable biomimetic technologies which respond quickly to stimuli.
For these reasons, this study focuses on using DEs for pattern generation in artificial cephalopod skins, and describes how they can be designed, mathematically modelled and analysed. We develop a hyperelastic electro-mechanical model and use it to simulate a variety of dynamic patterns using artificial chromatophores, composed of a linear array (thread) of DE. The DE thread is divided into cells, each of which is capable of sensing local or remote strain through strain-response switches, such as the DE switches developed by O'Brien et al. [9]. These sensors can be considered analogous to the strain-sensing spindle cells found in biological muscles [4]. By imposing simple rules that dictate the local interaction of neighbouring artificial chromatophores, we generate a range of complex and controllable patterns such as unidirectional propagation and sawtooth-like oscillation. We demonstrate that the generalized system design allows mimicry of dynamic skin patterning seen in biological chromatophore systems including the passing cloud display and other complex behaviour. This paper continues by outlining a generalized pattern generation system made from artificial chromatophore cells in §2. Section 3 then integrates this system into a time-dependent hyperelastic model that captures the physical and mechanical properties of DEs. This section also presents a discretization scheme of the model into linear arrays of artificial chromatophores or threads. Section 4 proceeds with the implementation of rule bases, using the discretization scheme, to generate a selection of biomimetic patterns including sawtooth oscillation, unidirectional propagation and complex behaviour. The paper concludes with a discussion of the potential of our system for future pattern generation in artificial skin.

Artificial chromatophore system design
The proposed artificial chromatophore system consists of a sheet of DE with a number of coloured electrodes painted on the material, forming a set of actuators or artificial chromatophore cells. The aim of our approach is to examine how patterns propagate within a chromatophore skin-like material, without the complexity of top-down independent control of each cell. We implement two types of cell, which differ in the control of their actuation. To reflect the intrinsic behaviour of chromatophores, self-sensing cells use control that is dependent on deformations within the material. Extrinsic behaviour is included using manual cells, which may control their actuation according to any external stimuli in the material environment (e.g. light or heat), or, indeed, external neural input.
As an illustrative example, consider a manual cell followed by a series of self-sensing cells as shown schematically in figure 2. If the manual cell is stimulated the deformation would be detected by a neighbouring selfsensing cell, prompting it to actuate, given the right control rule. This could potentially cause a 'chain-reaction' of actuation along the series of self-sensing cells, where the final state of the system has every cell actuated. Similarly, the contraction of the manual cell could also cause contraction neuro-electric signals radial muscle fibres sacs propagations that cause the chromatophore system to approach its original state. This simple example illustrates a mechanism of transitioning between two pattern types: 'all unactuated' to 'all actuated'.
We now assume that each cell is charged using an individual RC (resistor -capacitor) circuit connected in series to a source voltage. This models the finite resistance of the conducting electrodes of the DE actuator, coupled to its capacitance, as it charges and discharges during actuation. We equip each circuit with a binary switch, which immediately begins to charge the cell when closed. Self-sensing and manual cells differ in the control of this switch.
Self-sensing cells can be thought of as agents that decide when to contract and expand by choosing when to open and close their switch. To make such decisions the agents have two (local) pieces of information available to them: (i) the current state of their switch and (ii) their source, a scalar function of the current deformation. For example, the source may be the mean strain of the agent's cell, or the mean strain of its neighbouring cells. The agent closes its switch when the source lies within pre-defined ranges, forming an activation rule for the cell to expand. If the switch is closed and the source lies within a second set of ranges, the agent will open its switch, contracting the cell and forming the deactivation rule. Both rules together form the rule base of the system.
In practice, rule bases could be implemented using the dielectric elastomer switches (DES) proposed in [9]. By embedding DES close to a cell, we can use changes in resistivity to infer the cell's source and control its actuation through digital circuits. That being noted, an implementation of such a system is outside the scope of this paper, and for simplicity we will assume that perfect knowledge of the cell's strain is available for use in control.
The manual cells have much flexibility in their implementation and their control. In this paper, we use their actuation as triggers that prompt self-sensing cell switching to change the material pattern type. The cause of their actuation may be controlled by, for example, a toggle switch, or a photodiode in conjunction with a voltage gate.

Mathematical model
This section describes the operation of DE actuators and framework required to mathematically model the proposed DE pattern generation system. As introducing sets of switches increases system complexity, we discuss how we can characterize the state of a system using both the states of the switches and the motion of the material. We then introduce our viscoelastic continuum model, developed specifically for these DEs. This section concludes with a discussion of the model parameters fitted for 3M's VHB4905 [10], a highperformance, and commonly used, DE that we consider throughout the rest of the paper.

Dielectric elastomers as artificial muscles
When coated with a compliant electrode, sheets of DE store charge, much like a parallel plate capacitor. The build-up of charge induces Maxwell stresses within the material. For incompressible elastomers, this causes the material to expand in the plane of the electrodes and contract in the direction normal to the electrodes. For simplicity, we consider a situation in which the DE is constrained in one direction by light, stiff, parallel fibres added to the surface of the membrane in the plane of the electrodes [11]. The effect of the fibre constraint is to fix the length of the membrane in the direction of the fibres, so the expansion of the DE is uniaxial, normal to the fibre constraints, as illustrated in figure 3. It is known that a membrane subjected to uniaxial force along its length can achieve large deformation when the width direction is constrained [12]. Material undergoing such electrostatic expansion is defined as active and is passive otherwise. It has been shown [13] that the Maxwell stress, P, relates to the thickness of the elastomer, d, with applied voltage, V, by where 1 d ¼ 1 r 1 0 , 1 0 is the permittivity of free space (approx. 8.85 Â 10 212 Fm 21 ) and 1 r is the relative permittivity of the dielectric (typically around 4). We define active material that has approached steady state to be actuated and is unactuated otherwise.

Global and local dynamics
Including switches and sensors in the DE means that the motion of the material may be described by a dynamical system that is piecewise smooth [14]. Discontinuities occur at the crossing of thresholds, where the closing and opening of switches result in fast changes in system variables such as applied stress. To assess such discontinuities, we introduce the global state of the system, defined as the set of on/off states of all the switches. The system transitions between these discrete states when the material crosses a threshold. We define g i to be the state of switch i with g i [ f0, 1g, where 0 denotes the switch is open and 1 denotes the switch is closed. The global state of the system with n switches, one per cell, G is then defined to be G ¼ [g 1 , g 2 , . . ., g n ]. The global state does not uniquely define the material configuration, but is instead used to determine the electro-mechanical dynamics throughout the material.
The transition to the next global state, if any, is determined by the dynamics that occur when in a global state, described by a smooth dynamical system unique to that global state. We refer to this as the local dynamics. The initial conditions and governing dynamics of the global state determine which global state the system will transition to, as well as providing the initial conditions for the next global state.

Viscoelastic model
Our proposed viscoelastic model is based on that in [15] and has been chosen for its tractability and consistency with experimental data at high strains, for both fast and slow strain rates. Other viscoelastic models, such as in [16], are described using delay-differential equations (DDEs) and, although accurate, require sophisticated numerical integration schemes as well as being computationally expensive. Simpler models, such as [17], incorporate hyperelasticity into a Maxwell viscoelastic model [18]. Although these models are more tractable, studies indicated they were inaccurate for stretch ratios larger than 3 [19].
The DE is modelled using the hyperelastic Gent strain energy form [20]. The strain energy function, W, has the general form where m and J are material constants and I 1 is the first invariant of the deformation gradient written in terms of principal stretches in coordinate direction k, l k , as for an incompressible material. Because the DE we consider in this paper is fibre constrained in the width direction, 2, the principal stretch l 2 is a constant that we can choose (within limits), the prestretch l 2,pre . Hence we need only determine the behaviour of the material in the longitudinal, 1, direction. In practice, to prevent buckling and improve actuator performance, the DE is also prestretched in the longitudinal direction [21].
The Cauchy stress in the longitudinal direction, 1, of a generic hyperelastic element is given by where P is the Maxwell stress. We model the continuum as two parallel superposed networks, as in [15]. The first network, A, consists of a hyperelastic spring, while the second, B, consists of a viscous dashpot and a hyperelastic spring connected in series. Such rheology, summarized in figure 4, has been shown to successfully model time-dependent behaviour for DEs in numerous applications [22,23]. Both hyperelastic springs are modelled using the Gent strain energy form [20].
Following this rheology, we let j k and j e k represent the principal stretches in direction k over the viscous and elastic components, respectively, of network B, and l k the corresponding total stretch in the network (note, again, that the stretches in direction 2 are constant in our implementation). Using multiplicative decomposition, the elastic stretch is , meaning that larger viscous stretches result in lower elastic stretches.
The total stress is found by summing the contribution from the elastic components in each network. Let m A , m B , J A and J B represent the material constants in the hyperelastic components of networks A and B, respectively. Using (3.2), the Cauchy stress in direction 1 of the material is The first and second terms of (3.3) represent the stresses from network A and B, respectively. From (3.3) we note that, in general, larger viscous strains result in lower Cauchy stresses. Therefore, larger differences between l k and j k , which can arise from fast straining in the material, often result in greater resistances to motion. The rate of change of the viscous component in network B is found by modelling the dashpot in network B as a Newtonian fluid [15], yielding where h is the viscosity of the dashpot. For larger h, the rate of plastic deformation is slower, resulting in a longer relaxation time. The viscosity and spring modulus define the relaxation time, t ¼ h/m B [15]. In order to model the longitudinal prestretch, we consider the material to have been fixed to the desired ratio in the 1-direction, l 1,pre , before waiting a long time for the stresses in the viscous component to approach steady state. At steady state, we have j 1 ¼ 0, which implies that l 1 ¼ j 1 . The initial conditions for a prestretched material starting from rest are then defined as l 1 j t¼0 ¼ j 1 j t¼0 ¼ l 1,pre . Similarly, the fibre constraints imply that l 2 ¼ j 2 ¼ l 2,pre for all t.
Time-dependent behaviour is modelled using the constitutive equations (3.3) and (3.4) and requires the parameters  We are particularly interested in using our model to assess VHB4905, which has been shown to be capable of exhibiting linear actuation strains of over 500% [24], in the same order of magnitude as cephalopods' chromatophores [3]. Consistent with experimentation, this study uses proposed limits of material strength failure occurring when l ¼ 6 [10]. VHB4905 is manufactured with a thickness of 0.5 mm, which is a convenient thickness for designing artificial cephalopod skin. Since VHB4905 has a Poisson ratio of 0.49 [10], this study assumes it is fully incompressible. The material model parameters, given in table 1, have been fitted to existing experimental data [23].

Discretization of the viscoelastic model
In order to simulate and analyse our model, we discretize the material geometry. For simplicity, we consider the unidirectional motion of a membrane that is fibre constrained in the width direction. The model takes the form of a piecewise smooth set of ODEs, allowing for an efficient analysis of the dynamics.
The DE is discretized using n three-dimensional rectangular cuboids, each with mass m, which remain rectangular during deformation; this is illustrated schematically in figure 5. We consider elements with natural length L in the planar directions, and T in the thickness direction. Elements are joined end-to-end in series to form a thread of length l, with fixed boundary conditions, prestretched in the length direction by a factor l 1,pre , and in the width direction by the fibre constraints. Elements are represented by two nodes constrained to displace along the thread length only.
As the configuration of each element is determined by the elastic and viscous stretches in the length direction, defined as l 1 and j 1 in the previous section, we now relabel with subscripts denoting the element number. Thus, element i has stretch ratio l i with viscous stretch j i , both in the length direction. Electroded elements, with potential V i , are connected to individual RC circuits each with source voltage V S and resistance R. We define S j to represent the set of elements that are members of actuator j. Applying basic laws and geometry, V i evolves according to where and G(i) represents the global state The stress s i of an element in the thread length direction is calculated using (3.3) The force on node i, F i , is calculated using the stresses along the thread length in the adjacent elements i and i þ 1 yielding The evolution of the stretch ratios is then defined by the ODEs l i :: The rate of change of viscous stretch of element i is calculated using (3.4) as : Equations (3.5)-(3.7) form the equations of motion that define the system states V i , j i and l i , together with the boundary conditions which constrain the first and last nodes to be fixed.
An efficient way to solve for steady states is to assume activated electroded material to have a uniform stretch ratio, l A , implying all passive material will have stretch ratio l B . Let s elec l,V ð Þ and s mat l,j ð Þ represent the electrical stress and material stress, respectively, of a section of material with stretch l and viscous stretch j in direction 1 and potential V. If there are a active blocks and b ¼ l/L 2 a passive blocks, the fixed boundary conditions give an equation for l B , 0 ¼ (s elec (l A , V S ) À s mat (l A , l A ))l À1 A þ s mat (l B , l B )l À1 B , (3:8)

Pattern creation
We now use the modelling framework developed above to design and analyse controllable systems that are able to switch between different pattern types. We optimize parameters to accentuate the differences in pattern types, maximizing the change in area covered by cells in the actuated and unactuated configurations.
We analyse three pattern types, designed to mimic those seen in biological systems, summarized in table 2. For example, the type II pattern is capable of mimicking the passing cloud display of the cuttlefish, by allowing bands of actuation to propagate along the thread. Additionally, the type III pattern is designed to mimic signalling patterns seen in cephalopods, by choosing rules that result in complex activation sequences to propagate along the thread. Each pattern type is defined with a unique source and rule base.

Optimization of actuator size and spacing
Considering applications such as camouflage and signalling, it is beneficial to maximize the capability to display differences in perceived colour. This is achieved by considering the parameters that maximize the change in length covered by M cells, between when all cells are in actuated and unactuated states, illustrated in figure 6 for a particular combination of electroded and non-electroded cells.
In the length direction of the actuated configuration, let cells with stretch ratio l C and natural length d C be separated by sections of passive membrane with stretch ratio l P and natural length d P . The relative change in length covered by the cells is thus given by a ¼ (l 1,pre À l C )(l P À l 1,pre ) (l C À l P )l 1,pre , hence the optimal prestretch (which maximizes a) is given by and inter-element spacing Setting l C and l P to maximize the change is dangerous, since this may put the material close to material strength failure. It is therefore beneficial to leave linear strain margins for cells to actuate into. The buckling margin need not be as high as the material strength failure margin since the steady-state stretch in passive sections is minimized when all cells are fully actuated.
This study considers fabricating active cells that are 50 mm long at prestretch, as this has already been shown to be a plausible order of fabrication size when testing [3], although the system is scalable to other orders of magnitude. The maximal change in area that can be achieved is 42.02%, occurring when l C ¼ 6 and l P ¼ 1. However, to allow a margin against failure modes, the values l C ¼ 5 and l P ¼ 1.25 are used in this paper, which provide an area change of 33.33%. To assess whether this change in area is significant for a particular application, one would need to perform more detailed tests on patterns produced by a real membrane.
The voltage needed to drive the actuator is calculated using (3.8), which gives V s ¼ 5600 V. This maximizes the localized steady-state area change when all cells are actuated, as fewer actuated cells reduce the stress in the passive membrane. Finally, we set the RC circuits to have a resistance of 10 MV, which results in fairly fast charge and discharge times.
The full list of parameters used throughout the rest of this paper is given in table 3. Note that the patterns generated are not unique to this set of parameters but would result in either reduced relative area change or lower failure tolerances.  The type I rule base and source was configured to enable cells to enter a sawtooth-like oscillation, as shown in figure 7, a shape that reflects the hyperelastic stress response of DE. Such patterns could be useful for signalling applications. We illustrate the dynamics associated with this rule base using the simplest combination of cells: a single selfsensing cell, in the middle of a longer ( passive) thread. To ensure oscillatory behaviour, limits on r on and r off are chosen to be within the steady-state source of the cell in global states G ¼ 1 and G ¼ 0, respectively, requiring l l,pre , r on , r off , l C . Greater differences in the thresholds increase the amplitude of oscillation, as the source oscillates (approximately) between them.
Amplitude parameters of r on ¼ 2.9, r off ¼ 4.2 are studied, for which the amplitude of oscillation is approximately l 1,pre . Applying source voltage to the cell at t ¼ 0 (shown in figure 7) clearly triggers a limit cycle that results in significant optical modulation, with a period of approximately 3.7 s. The sawtooth-like shape occurs because the cell takes longer to expand than contract, a consequence of working against material stresses which arise from prestretch. Similarly, these material stresses increase cell contraction speed, causing the frequency of oscillation to be dependent on the amplitude.
During oscillation, the momentum of the material does not significantly force the material beyond the thresholds, implying viscous forces and electrical gradients act quickly to change the deformation direction.

Type II pattern: unidirectional propagation
When a number of self-sensing cells are placed after a single manual cell, the type II source is capable of propagating sequences of activation along the thread length. By inputting activation sequences into the manual cell, patterns such as the 'passing cloud' display of the cuttlefish can be mimicked. To achieve this, the corresponding rule base was designed to sense the state of the neighbouring cell.
The limit on the activation strain r on needs to be large enough to avoid activating cells in the initial steady state G ¼ 0. This limit should also be small enough to detect the activation of the neighbouring cell in the extreme circumstance when the steady-state source is minimized in state G ¼ 1. Together, these imply l 1,pre , r on , l C . Similarly, the limit for r off must be larger than the maximal steady-state source when no cells are actuated, to ensure deactivation is detected. Additionally, r off is set below r on to avoid local oscillations. Therefore, the necessary rule base conditions for propagation are l 1,pre , r off , r on , l C : To illustrate propagation, we consider a system with one manual cell (cell 1), followed by 29 self-sensing cells (cells 2-30), and a near-optimal inter-element spacing of natural length d 2 ¼ 40 mm. As cell contraction is faster than actuation, the rule base parameters r on ¼ 3.7 and r off ¼ 2.6 were chosen because of their fairly even activation and deactivation propagation speeds. The manual cell is set to activate and deactivate every 1 s, ceasing after 8 s, to create regular activation and deactivation propagations. Figure 8 illustrates regular propagations of activation along the full length of the thread. As tension in the thread causes activation propagations to travel slower than deactivation propagations, cells further away from the manual cell are activated for a shorter length of time. Therefore, the activation time of the manual cell effectively controls the propagation distance of the wave. The approximately regular activation cycles imply that viscous effects do not significantly affect the travelling wave.  Table 3. System parameters used in the model. The change can be seen more clearly in the global state transition diagram (figure 8b), clearly illustrating the width of the travelling actuation waves. The pattern in figure 8 was designed to mimic the passing cloud display [4], where the propagation of actuation is analogous to the bands of colour travelling across the body of a cuttlefish, as shown in figure 8c.

Type III pattern: complex behaviour
Our artificial chromatophore skin system is also capable of creating complex, chaotic and long-lasting sequences of actuation, analogous to those produced in some biological systems [4]. We illustrate this by using the type III source. Here, the time evolution of a cell is based on the number of actuated neighbours. If none or two neighbouring cells are actuated, the agent deactivates its cell. If one neighbour is actuated, the agent activates its cell. To implement such behaviour, the corresponding rule base parameters are set to infer the number of actuated neighbours: r off,1 and r off,2 detect zero or two actuated neighbours, respectively, and r on,1 and r on,2 detects one actuated neighbour. The thresholds should satisfy r off,1 , r on,1 , r on,2 , r off,2 , in order to partition ranges for detecting the number of actuated neighbours.
The limits for these rule base parameters are thus determined by analysing the steady-state source as a function of the number of activated cells, using (3.8), for zero, one and two locally activated neighbours. If the ranges of the steadystate source for each neighbour are sufficiently distinct, rule base limits are imposed on the extremes of the ranges to ensure a response to the number of activated neighbours. However, a cell may not respond if its neighbour changes its activation state before the cell approaches the threshold.
An example of a type III pattern is illustrated in figure 9. Here, the manual switch (cell 1) is used to trigger patterns by activating for 0 s t 1 s. The pattern creates complicated and long-lasting switching behaviours which occur because the switching rule makes it difficult for the thread to reach steady state. In the first 0.6 s, an initial activation propagation similar to type II behaviour is observed as cells respond to an activated neighbour. For t . 0.6 s complex behaviour emerges.

Conclusion
In this paper, we have presented the design of biomimetic dynamic pattern generation systems in soft artificial skin made from DE. By considering linear arrays of DE actuators, we have demonstrated that such systems can exhibit a range of patterned behaviour using simple local rules to control actuation. The basis of our analysis was a viscoelastic continuum model that incorporates the hyper-viscoelastic properties of DEs, optimized to accentuate actuation. Our results suggest that, although pattern generation is not always regular, the system is capable of generating multiple dynamic pattern types analogous to those seen in biological systems, with an expansion ratio of the same order of magnitude as real cephalopods [4]. Patterns include simple sawtooth-like oscillation, 0   30  29  28  27  26  25  24  23  22  21  20  19  18  17  16  15  14  13  12  11  10  9  8  7  6  5  4  3  unidirectional propagation of actuation sequences mimicking the passing cloud displays in cephalopods, and more complex long-lasting behaviours. The artificial cephalopod skin system presented here may be extended to applications such as camouflage and signalling. Following the suggestion that some cephalopods communicate using polarized light [25], patterns could be further accentuated by embedding layers of polarizing filters onto the surface of the DE. By strategically varying the direction of the filter along the surface of the thread and layering multiple actuator membranes on top of each other, we could use the dynamics generated here to create complex patterns of light. This would be particularly useful for search and rescue applications, where rescuers need to stand out. To assert the compatibility of pattern propagation speeds for a particular application, further investigation into the accuracy of actuation speeds predicted by the model is required.
Future work will consider adjusting system parameters to improve propagation control and generate new patterns using other local rules, conducting a more extensive analysis of the different pattern types that can be obtained under variation of system parameters, as well as extending the model to simulate patterns in two-dimensional array systems. This is expected to yield further varieties of dynamic patterns, some of which may correspond to those in the natural world.