Finite Element Modelling of a Cellular Electric Microenvironment

Summary

This paper presents a strategy for building finite element models of fibrous conductive materials exposed to an electric field (EF). The models can be used to estimate the electrical input that cells seeded in such materials receive and assess the impact of changing the scaffold's constituent material properties, structure or orientation.

Abstract

Clinical studies show electrical stimulation (ES) to be a potential therapy for the healing and regeneration of various tissues. Understanding the mechanisms of cell response when exposed to electrical fields can therefore guide the optimization of clinical applications. In vitro experiments aim to help uncover those, offering the advantage of wider input and output ranges that can be ethically and effectively assessed. However, the advancements in in vitro experiments are difficult to reproduce directly in clinical settings. Mainly, that is because the ES devices used in vitro differ significantly from the ones suitable for patient use, and the path from the electrodes to the targeted cells is different. Translating the in vitro results into in vivo procedures is therefore not straightforward. We emphasize that the cellular microenvironment's structure and physical properties play a determining role in the actual experimental testing conditions and suggest that measures of charge distribution can be used to bridge the gap between in vitro and in vivo. Considering this, we show how in silico finite element modelling (FEM) can be used to describe the cellular microenvironment and the changes generated by electric field (EF) exposure. We highlight how the EF couples with geometric structure to determine charge distribution. We then show the impact of time dependent inputs on charge movement. Finally, we demonstrate the relevance of our new in silico model methodology using two case studies: (i) in vitro fibrous Poly(3,4-ethylenedioxythiophene) poly(styrenesulfonate) (PEDOT-PSS) scaffolds and (ii) in vivo collagen in extracellular matrix (ECM).

Introduction

ES is the use of EFs with the aim of controlling biological cells and tissues. Its mechanism is based on the physical stimulus transduced to the cell when the biomolecules inside and surrounding it are exposed to an externally generated voltage gradient. Charged particles are engaged in an organized motion governed by Coulomb's law, generating drag forces upon the uncharged particles. The resulting fluid flow and charge distribution alter cell activities and functions such as adhesion, contraction, migration, orientation, differentiation and proliferation¹ as the cell attempts to adapt to the change in the microenvironmental conditions.

As EFs are controllable, non-invasive, non-pharmacological and shown to have an effective impact on essential cell behavior, ES is a valuable tool for tissue engineering and regenerative medicine. It has been successfully used to guide neural², skeletal³, cardiac muscle⁴, bone⁵ and skin⁶ development. Moreover, as it enhances iontophoresis⁷, it is used as an alternative or complementary treatment to conventional pharmacological ones. Its efficiency in pain management is still debated as higher quality clinical trials are awaited^8,9,10. Nevertheless, no adverse effects were reported and it has the potential to improve patient welfare^{11,12,13,14,15}.

While only clinical trials can give a definitive verdict for the efficacy of a procedure, in vitro and in silico models are required to inform the design of predictable ES treatment as they offer stronger control over a wider range of experimental conditions. The investigated clinical uses of ES are bone regeneration^16,17, recovery of denervated muscles^18,19, axonal regeneration after surgery^20,21, pain relief²², wound healing^23,24,25 and iontophoretic drug delivery²⁶. For ES devices to be widely introduced on all possible target applications, clinical trials have yet to establish stronger evidence for efficient treatment. Even in domains where both in vivo animal and human studies consistently report positive outcomes, the great number of reported methods coupled with too little guidance on how to choose between them and high acquisition price deters clinicians from investing in ES devices²⁷. To overcome this, the target tissue can no longer be treated as a black box (limit of in vivo experiments) but must be seen as a complex synergy of multiple subsystems (Figure 1).

Multiple ES experiments have been carried out in vitro over the years^{28,29,30,31,32,33,34}. Most of these only characterize the ES through the voltage drop between the electrodes divided by the distance between them - a rough approximation of the electric field magnitude. However, the electric field itself only influences charged particles, not cells directly. Also, when multiple materials are interposed between the device and the cells, the rough approximation may not hold.

A better characterization of the input signal requires a clear view on how the stimulus is transduced to the cell. Main methods of delivering ES are direct, capacitive and inductive coupling^35,36. Devices for each method differ with electrode type (rod, planar or winding) and placement relative to the target tissue (in contact or isolated)³⁵. Devices used in vivo for longer treatments need to be wearable, thus the electrodes and most times the energy source are either implanted or attached to the skin as wound dressings or electroactive patches. The generated voltage gradient displaces charged particles in the treatment area.

As it impacts the resulting charged particle flow in the vicinity of the cells, scaffold structure is of utmost importance in the design of ES protocols. Different charge transport configurations arise if the platform material, synthesis technique, structure or orientation relative to the voltage gradient change. In vivo, the availability and movement of charged particles is impacted not only by cells but also by the collagen network and interstitial fluid composing the supporting ECM. Engineered scaffolds are increasingly used to better recreate natural cell microenvironments in vitro^1,35. Concurrently, the ECM is a complex natural scaffold.

Artificial scaffolds are based on metals, conducting polymers and carbon, engineered with a focus on balancing biocompatibility with electrochemical performance and long-term stability³⁶. One versatile scaffold type is the electrospun fibrous mat that offers a controllable nanoscale topography. This can be engineered to resemble the ECM, thus deliver similar mechanical cues that aid regeneration of a wide range of tissues³⁷. To significantly impact ES, the mats need to be conductive to some degree. However, conductive polymers are difficult to electrospin and blending with insulating carriers limits the conductivity of the resulting fibers³⁸. One solution is polymerizing a conductive monomer on the surface of a dielectric fiber, resulting in good mechanical strength and electrical properties of the end product³⁸. An example is coating silk electrospun fibers with the semi conductive PEDOT-PSS³⁹. The combination of mechanical and electromagnetic cues significantly accelerates neurite growth^40,41,42. Neurites follow scaffolds fibers alignment, and elongate more after exposure to an EF parallel to the fibers than to a vertical one⁴³. Similarly, alignment of fibrous scaffolds to the EF also promotes myogenic maturation³³.

The ECM is mainly composed of fibrous-forming proteins⁴⁴, out of those collagen type I being the major constituent in all animal tissues apart from cartilage (rich in collagen type II)⁴⁴. Tropocollagen (TC), triple helical conformation of polypeptide strands, is the structural motif of collagen fibrils⁴⁵. Transmission electron microscopy and atomic force microscopy images of collagen fibrils show a D-periodic banded pattern⁴⁶ explained by the Hodge & Petruska model⁴⁷ as regular arrays of TC gaps and overlaps⁴⁵. Tendons are composed of an aligned collagenous fibrillar matrix shielded by a non-collagenous highly hydrophilic proteoglycan matrix^48,49. Decorin is a small leucine-rich proteoglycan (SLRP) able to bind the gap regions of collagen fibrils and connect with other SLRPs through their glycosaminoglycan (GAG) side chains⁴⁹. Studies done on tendons show that their electrical properties change significantly when hydrated^50,51, charge transport mechanism changing from protonic to ionic as hydration level increases⁵¹. This suggests that electric conduction along a collagen type I fibril could be enabled by a Decorin-water coat, with gap and overlap regions having different electrical conductivities and dielectric constants.

As identical recreation of the ECM by artificial scaffolds is improbable, the knowledge producing synergy between in vivo and in vitro enabled by translatable results seems to be at a dead end. In silico modelling not only re-enables translation between the two, but also adds important benefits in characterizing the unknown processes involved in ES. Comparing the in vivo observations with the in vitro can bring information on the coupling strength between the target tissue and the rest of the organism but does not uncover current knowledge limits. The unknown can be exposed by observing the difference between what is expected to happen based on the current knowledge and what happens. In silico experiments based on mathematical modelling allow splitting the process into known and unknown subprocesses. This way, phenomena not accounted for in the model come to light when in silico predictions are compared to in vitro and in vivo experiments.

Forming and testing hypotheses regarding the underlying mechanism(s) of how cells and tissues are affected by electrical fields is hindered by the great number of parameters⁵² that need to be tested separately. To define representative experimental conditions, the ES process must be split in subprocesses (Figure 1) and dominant input signals affecting cell behavior must be identified. Models representing fundamental physical effects of ES on cells describe the domain that couples the EF with the cell - that of charged particles⁵³. The behavior of particles exterior to the cell depends on the microenvironment and can be investigated separately from the cell. The dominant input signal for the cell is the subset of ES device outputs that causes the greatest degree of variability in the cell response. The smallest subset of the full experimental parameters that can generate variations in all the dominant cell input signals can be used to decrease the parameter space dimension and the number of test cases.

The input of the biological ES target model must be a subset of the output signals produced by the ES device that are useful in describing the physical effects of ES on cells. A simple bioreactor with direct coupling has the same structure as electrolytic electrochemical cells. Models of those show the primary (accounting for solution resistance), secondary (also accounting for faradic reactions) or tertiary (also accounting for ion diffusion) current density distribution. As complexity translates into computational cost, the simplest model is most suitable for parameter space explorations. Simulations of fibrous composites motivated by material properties⁵⁴ focus on bulk material properties as a result of complex micro-architecture, hence cannot describe local effects of EF exposure. Existing in silico models, motivated by ES, focus on the biological sample, be it a single cell immersed in a homogenous medium^55,56,57, or complex tissues with homogenous extracellular space⁵⁸. Charge and current density (Figure 2) can act as interface signals between models of the ES device and the biological sample, or between different components of the ES device. The proposed FEM based protocol uses the equations described in Figure 2 and was used to study how scaffold dependent parameters can be used to modulate those two signals, independent of the EF generated by a direct coupling setup. Results stress that it is necessary to account for scaffold or ECM electrical properties when investigating how ES impacts target cells.

Protocol

1. Build the model in COMSOL

1. Open COMSOL and select Blank Model.
2. Parameters: In Model Builder, right click on Global Definitions, select Parameters, and add parameters according to Table 1.
3. Materials: Add materials with properties according to Table 2.
   1. In the Model Builder under Global Definitions, right click Material and select Blank Material.
   2. Add Material Properties: In the settings of the newly added material, expand Material Properties > Basic Properties, select Relative Permittivity and press the '+' symbol to add property. Repeat the step for Electrical Conductivity.
   3. In the Material Contents, fill in the current material properties according to Table 2.
4. Component definition: In the Home tab of the top ribbon, left click Add Component and select 3D. A new component node will appear in the Model Builder.
5. Geometry: In the Model Builder, right click Geometry, left click insert, double click on the Full Model and select the appropriate sequence (SC/SNC/RC/RNC/RNCd).
6. Alternatively, build the geometry from scratch following the subordinate steps.
   1. In the Model Builder, left click Geometry under the current component node. The same Geometry node will be referred in the following subordinate steps. In the Settings window, change the Length unit to nm and make sure "Scale values when changing units" is left unchecked.
   2. Geometry of scaffold
      1. In the Model Builder, right click on Geometry and select Block.
      2. In the Settings tab, fill in the Label "Scaffold extents".
      3. Expand Size and Shape and fill in "Ws", "Ls", "Hs" for Width, Depth and Height.
      4. Expand Position and change the Base to Center.
      5. Expand Selections of Resulting Entities, select New and add "Scaffold Selection".
   3. Extents of surrounding substance: In the Model Builder, under Geometry, right click Scaffold extents and select Duplicate, and edit Settings tab.
      1. Change the label to "Media extents".
      2. To each box in the Size and Shape section, add "med_ratio *" before the existent parameter.
      3. In the Selections of Resulting Entities, add New selection as "Media Selection".
   4. Geometry of a fiber (SC skip to 1.6.6)
      1. Core (RC skip to 1.6.4.2)
         1. In the Model Builder, right click Geometry, select Cylinder, and edit Settings tab as follows.
         2. Change the label to "Core".
         3. Expand Size and Shape and fill in "Rc" and "Lf" for Radius and Height.
         4. Expand Position and fill in "-excess*Ws/2" and "-Lf/2" for x and y.
         5. Expand Axis and change Axis type to y-axis.
         6. Expand Selections of Resulting Entities and add "Core Selection".
      2. Coat (SNC skip to 1.6.5)
         1. Slab (RC or RNC skip to 1.6.4.2.2; RNCd skip to 1.6.4.2.3)
            1. In the Model Builder under Geometry, right click Scaffold extents and select Duplicate. Make sure the duplicate (Scaffold extents 1) comes right after Core in the geometry sequence.
            2. In the Settings tab, change label to "Coat".
            3. In the Settings tab, expand Selections of Resulting Entities and create "Fiber Selection".
            4. Skip to 1.6.5.
         2. Homogenous
            1. In the Model Builder right click Geometry, select Cylinder and edit its Settings tab as follows.
            2. Change the label to "Coat".
            3. Expand Size and Shape and fill in "Rf" and "Lf" for Radius and Height.
            4. Expand Position and fill in "-excess*Ws/2" and "-Lf/2" for x and y.
            5. Expand Axis and change axis type to y-axis (RC skip to 1.6.4.2.2.7).
            6. Expand Layers, fill in "Rf-Rc" for Thickness and leave only "Layers on side" checked.
            7. Expand Selections of Resulting Entities and add "Coat Selection" or, for RC, "Fiber Selection" (RC skip to 1.6.5).
            8. In the Model Builder, right click Geometry and select Delete Entities.
            9. In the Settings tab, change Geometric Entry Level to Domain and choose Coat Selection for Selection.
            10. Remove all domains but domain 3 from the selection box (the core of the layered cylinder).
            11. Expand Selections of Resulting Entities and select Coat Selection.
            12. Skip to 1.6.5.
         3. Periodic array of two types
            1. In the Model Builder, right click Geometry and select Cylinder. Right click the newly created element and select Duplicate. Edit their Settings tab as follows.
               1. Set Label to "Coat 1" / "Coat 2".
               2. Set radius to Rf.
               3. Set height to "D*prop"/"D*(1-prop)".
               4. Expand position and fill in "-excess*Ws/2" for both on x and "-Lf/2"/"-Lf/2+D*prop" on y.
               5. Expand Axis and change axis type to y-axis.
               6. Expand Layers section and fill in "Rf-Rc"; Make sure only the "Layers on side" box is checked.
               7. Expand Selections of resulting entities and create "Coat 1 Selection"/ "Coat 2 Selection".
            2. In the Model Builder, right click Geometry and select Delete Entities.
               1. In the Settings tab, change Geometric Entry Level to Domain and select cylinder 2 (cyl2) and cylinder 3 (cyl3) from the graphics window.
               2. Remove all domains but domain 3 from the selection box (the core of the layered cylinders).
               3. Expand Selections of Resulting Entities and create Coat Selection.
            3. In the Model Builder, right click Geometry, expand Transforms and select Array. Edit Settings tab as follows.
               1. Select "Coat Selection" for Input objects.
               2. Set Array type to Linear.
               3. Set Size to "Lf/D".
               4. Set Displacement on y axis to "D".
   5. Geometry of an oriented fiber array
      1. Core Array (RC skip to 1.6.5.2)
         1. In the Model Builder, right click Geometry, Expand Transforms and select Rotate.
            1. In the Settings tab, turn on the Input objects Switch and select Core Selection from the dropdown list. Make sure Keep input objects box is not checked.
            2. Expand Rotation Angle and fill in "theta" for Rotation.
            3. Expand Selections of Resulting Entities and create "Fiber Selection".
         2. In the Model Builder right click Geometry, Expand Transforms and select Array.
            1. Change label to "Core Array".
            2. Select Core Selection in the dropdown list for Input objects.
            3. Expand Size, change Array type to Linear and fill in "n_1*(theta<=45) +n_2*(theta>45)" for Size.
            4. Expand Displacement and fill in "2*tes*Rc*cos(theta)", "2*tes*Rc*sin(theta)" for x and y.
            5. Expand Selections of Resulting Entities and select "Fiber Selection".
      2. Coat Array (SNC skip to 1.6.5.3)
         1. In the Model Builder, right click Geometry, Expand Transforms and select Rotate.
            1. In the Settings tab turn on the Input objects Switch and select Coat Selection from the dropdown list. Make sure "Keep input objects box" is not checked.
            2. Expand Rotation Angle and fill in "theta" for Rotation.
            3. Expand Selections of Resulting Entities and select "Fiber Selection".
         2. In the Model Builder, right click Geometry, Expand Transforms and select Array.
            1. Change label to "Coat Array".
            2. Select Coat Selection in the dropdown list for Input objects.
            3. Expand Size, change Array type to Linear and fill in "n_1*(theta<=45) +n_2*(theta>45)" for Size.
            4. Expand Displacement and fill in "2*tes*Rc*cos(theta)", "2*tes*Rc*sin(theta)" for x and y.
            5. Expand Selections of Resulting Entities and select "Fiber Selection".
      3. Cut Scaffold
         1. Unite Fibers: In the Model Builder right click Geometry, Expand Booleans and Partitions and select Union.
            1. Change label to "Fibers".
            2. Turn on Input objects switch and select Fiber Selection from the dropdown list.
            3. Make sure "Keep input objects" box is not checked.
            4. Make sure "Keep interior boundaries" is checked.
            5. Expand Selections of Resulting Entities and select Scaffold Selection from the dropdown list.
         2. Perform cut: In the Model Builder right click Geometry, expand Booleans and Partitions and select Intersection.
            1. Change label to "Scaffold".
            2. Turn on the Input objects switch and select Scaffold Selection from the dropdown list.
            3. Make sure "Keep input objects" box is not checked.
            4. Make sure "Keep interior boundaries" is checked.
   6. Ensemble Geometry
      1. Define surrounding substance geometry.
         1. In the Model Builder, right click Geometry, expand Booleans and Partitions and select Difference. Continue to the Settings tab.
         2. In the Objects to add dropdown list, select Media Selection.
         3. In the Objects to subtract dropdown list, select Scaffold Selection.
         4. Make sure "Keep input objects" and "Keep interior boundaries" boxes are both checked.
      2. Delete Media extents box.
         1. In the Model Builder, right click Geometry and select Delete entities. Continue to the Settings tab.
         2. Select Domain from the Geometric entry level dropdown list.
         3. Select Media Selection from the Selection dropdown list.
         4. Remove "dif1" from the selection box.
   7. Build geometry
      1. Under the Geometry node, make sure the last element before the Cumulative Selections node is Form Union and not Form Assembly. If necessary, right click Form Assembly and change Action to "Form a union" in the Settings tab.
      2. In the Model Builder, left click Geometry and select Build All.
7. Materials
   1. In the Model Builder under the Current component node, right click Materials and select Material Link.
   2. Repeat section five times for as many times as many materials are considered according to the complexity level.
   3. Associate materials for each component in the following order: surrounding substance, coats, cores.
   4. In the Settings tab, choose Media/Scaffold (for SC, SNC, RC)/Coat/Core Selection from the Selection dropdown list.
   5. Expand the Link settings and choose appropriate material from the dropdown list.
8. Meshing
   1. In the Model Builder, left click the Mesh node under the current component.
   2. In the Settings tab, select Normal from the Element size dropdown list and left click Build All.
9. Physics: In the Model Builder, left click current component, select Add Physics, expand the AC/DC module in the Add Physics tab, select the Electric Currents (ec) module and click Add to Component.
   1. Boundary conditions: In the Graphics tab, select the xy view.
      1. Ground: In the Model Builder, right click on the Electric Currents node and select Ground.
      2. Make sure the selection switch for the Boundary Selection is Active and left click on the highest surrounding substance face parallel to the xz plane to add boundary 5 in the Boundary Selection Box.
      3. Terminal: In the Model Builder, right click on the Electric Currents node and select Terminal.
      4. Make sure the selection switch for the Boundary Selection is Active and left click on the lowest surrounding substance face parallel to the xz plane to add boundary 2 in the Boundary Selection Box.
      5. Expand the Terminal section and select Voltage in the Terminal type dropdown list; fill in V0 for Voltage.
   2. Add complexity if needed: Save current component and go back to Step 1.4 to follow another complexity path (SC/SNC/RC/RNC/RNCd).

2. Perform simulation

1. Adaptive mesh refinement
   1. In Model Builder, left click the model root node and select Add Study to open the Add Study tab, select Stationary Study and right click Add Study button.
   2. Create a study step for each created component: under Study node right click Step1: Stationary Study and select Duplicate.
   3. For each study step modify the settings tab with the details of the corresponding component.
      1. Expand Physics and Variables Selection; in the Solve for column leave only the current component checked.
      2. Expand Study Extensions and check Adaptive mesh refinement box.
      3. Select current component's geometry from the dropdown list next to Adaptation in geometry.
   4. In Model Builder right click Study 1 and select Compute. This will generate adapted meshes for all the components' geometries with their current orientation angle.
2. Set orientation angle and perform a stationary study.
   1. In Model Builder, under Global definitions, left click Parameters and change parameter theta to the fiber orientation angle desired for simulation.
   2. For each component, in the Model Builder expand the component's node, right click Geometry and select Build all.
   3. In Model Builder, right click the adaptive mesh refinement study and select compute.
   4. In Model Builder left click the model root node and select Add Study to open the Add Study tab, select Stationary Study and right click Add Study button.
   5. In Model Builder, under the newly added study, left click Step 1, expand mesh selection and, for each component, select the meshes generated in the adaptive mesh refinement study.
   6. Right click Compute button.
3. Set orientation angle, input signal and perform time dependent study.
   1. In Model Builder, under Global definitions, left click Parameters and change parameter theta to the fiber orientation angle desired for simulation.
   2. In the Model Builder, left click the model root node and select Add Study to open the Add Study tab. Select Time Dependent Study, left click Add Study button and edit Settings tab as follows.
      1. Set Times to "range(0, (2*pi/omega)/39,2*pi/omega)".
      2. Expand Physics and Variables Selection; leave only the simulation component "Solve for" check box checked.
      3. Expand Mesh selection and select a mesh for the simulation component. Change the other components' meshes to No mesh.
      4. Expand Study Extensions and check the Adaptive mesh refinement box; select the simulation component's geometry from the dropdown list.
   3. Define time dependent boundary condition.
      1. Under the current component node right click Definitions and under Functions select Waveform.
      2. In the Settings tab change Function name to "Input".
      3. Expand Parameters, set Type to Sine, Angular frequency to "omega", Amplitude to V0.
      4. For the simulation component, under its Electric Currents node, select Terminal and change Voltage to "Input(t[1/s])".
   4. Right click Compute button.

3. Analysis

1. Charge density
   1. In the Model Builder right click Results node, select 3D Plot Group and edit Settings.
      1. Change label to "Charge density".
      2. Expand Data and select the Parametric study dataset in the Data set dropdown list.
      3. Expand Color Legend and check "Show legends" and "Show maximum and minimum values".
   2. Right click Charge density under the Results node in the Model Builder and select Volume; proceed to edit Settings tab.
      1. Expand Data and select "From parent" from Data set dropdown list.
      2. Expand Expression and fill in "ec.rhoq" in the Expression box.
      3. Expand Range and check Manual color range box.
         1. Set minimum to "-0.03" and maximum to "0.03".
      4. Expand Coloring and Style, set Coloring to Color table, Color table to Wave, check Color legend box, check Symmetrize color range.
      5. Right click Volume in Model Builder and select Filter.
         1. In the Settings tab fill in "abs(ec.rhoq)>0.012" in the Logical expression for inclusion.
   3. Left click on the Plot button to visualize results in the graphics window.
2. Current density
   1. In the Model Builder right click the Results node, select 3D Plot Group and edit the Settings tab.
      1. Change label to "Current density".
      2. Expand Data and select the Parametric study dataset in the Data set dropdown list.
      3. Expand Color Legend, check "Show legends" and "Show maximum and minimum values".
   2. Right click Current density under the Results node in Model Builder and select Arrow Volume; proceed to edit Settings tab.
      1. Expand Data and select "From parent" from Data set dropdown list.
      2. Expand Expression and fill in "ec.Jx", "ec.Jy", "ec.Jz" in the Expression boxes for the x, y and z components respectively.
      3. Expand Arrow positioning and fill in 20 for all coordinates number of points.
      4. Expand Coloring and Style, set Arrow length to Normalized, Arrow base to Center, check Scale factor and set it to 85.
      5. Right click Arrow Volume in Model Builder and select Color Expression.
         1. In the Settings tab fill in "ec.normJ" in the Expression box.
         2. Expand Coloring and Style, set Color table to Traffic, check Color legend and Reverse color table.
   3. Click on the Plot button to visualize results in the graphics window.

Results

The proposed model describes features of a composite mat with parallel fibers, immersed in a conductive substance and exposed to an externally generated electric potential gradient. Simulations show that accounting for the different components of a scaffold is important on a microscale and explore how change in alignment angle (input signal) of the fibers to the EF can generate variability in the current and charge density (output signals) in the vicinity of the fibers.

Five different geometrical complexity stages are presented, each having an effect upon the simulation result: smooth conductive slab (SC), smooth slab with non-conductive embedded fibers (SNC), rough conductive composite (RC), rough composite with non-conductive embedded fibers (RNC), rough composite with non-conductive embedded fibers and two types of periodic coating (RNCd) (Figure 3). Section 1.5 of the protocol presents the steps to importing the geometries in a project and section 1.6 shows how to build those step by step. The first two models do not account for surface morphology. SC and RC do not account for the fiber core dielectric properties. The RNC is the proposed model for nanofibrous artificial scaffolds, while RNCd is the proposed model for an ECM segment.

Minimization of computational cost was accomplished by reducing the ES device geometry to a model unit volume representing the microenvironment. While an ES device and scaffold's width and length can easily be at the order of a few centimeters, the containing fibers' diameter is usually lower than a micron. Here, we use a scaffold cut comparable to the fiber diameter to reduce the computational cost induced by the aspect ratio and highlight the effect of the scaffold's fibrous nature on the electric microenvironment. The rest of the ES device is replaced with electric potential boundary conditions chosen so that a rough approximation for the magnitude of the electric field is 100 V/m, a frequently reported stimulation parameter. Moreover, a unit volume with five parallel fibers - as the one used in simulations, presented in Figure 3 - is assumed to be representative of a whole planar fibrous mat. Three types of fibers can be distinguished in a 1D array: interior central (with the longitudinal symmetry plane of the scaffold splitting them in half), interior transitory (with lateral surface surrounded by other fibers but with asymmetrical sides), and exterior (at the edge of the scaffold). Five is the minimum number of fibers required in order to include all the three types defined.

The model mesh element size requires special attention as it may impact simulation results and thus fail to expose important effects (Figure 4). This is a general rule of the finite element method and an implication of the Nyquist-Shannon sampling theorem. The faster the essential simulation signals fluctuate in space the smaller the mesh elements need to be to produce a loyal representation of the phenomenon. On the other hand, the smaller the element, the greater the total number of model building blocks and the computational cost. The adaptive mesh refinement set up in section 2.1 is a good and facile method to balance those opposing objectives by decreasing the element size only where and as long as this operation produces a significant change.

A model that is too simplistic can fail presenting important effects (Figure 5,6). Simulations show that accounting for surface morphology and scaffold component electrical properties is not redundant in predicting electric microenvironments. While surface morphology has a direct impact on the stationary EF (compare SC and SNC with RC, RNC and RNCd), a comparison between RC and RNC predictions shows that nonconductive fiber cores amplify this effect. From the point of view of modelling cellular electric microenvironments on nanofibrous scaffolds, the SC, SNC and RC models are thus sub-optimal. However, it is good practice to incrementally add complexity as comparisons between the different stages help indicate what features give rise to specific effects.

Model complexity impacts current and charge density change with fiber alignment to the EF. The proposed protocol helps highlight the effect (Figure 5,6). While the SC model shows no variation in the proposed metrics when its alignment to the electric potential gradient is changed, the RNC model simulations predict a powerful contrast between the mat unit with fibers aligned to the EF and the one with fibers perpendicular (Figure 7). When the non-conductive cores come in the way of the current flow, they form periodic dams that lead to alternating regions of high and low charge density.

Dynamic ES regimes can be simulated with time dependent studies. Videos in supplementary files show predictions made for a sinusoidal input voltage on a full artificial scaffold model (RNC), with fibers parallel or perpendicular to the electric potential gradient. Small currents along the fibers perpendicular to the EF appear when charge is released from the scaffold as the EF magnitude decreases. This shows that stimulation could occur not only while the external EF is present, but also right after it is disconnected - See supplementary files for videos.

Figure 1: Hierarchical block diagram of modelling - advantages and limitations of modelling with in vivo and in silico models. Block color marks blocks on the same hierarchical level. Lower rank blocks are included in higher rank ones. Block stroke color marks possibility to include the block into a certain type of model - coupling with other system blocks do not have yellow in their stroke, as they are not components for in vitro models. Bullets act like valves and signify controllability of the block. When a valve is ON, signal can pass through all arrow paths in the subordinate subsystems that have the color of the valve in their stroke. Interpretation of the diagram: the ES process is composed of the stimulation device and biological target, each with several inter-connected deterministic or stochastic sub-processes that cannot be separated in vivo or in vitro, thus they have no red or yellow valve. Stochastic processes also intervene on the interface between the simulation device and biological sample when they are both stimulated. An in vitro model decouples the system of interest (i.e., skin segment) from the rest of the organism. Thus, only intrinsic processes of the system of interest topped by stochastic processes of different nature can be observed. However, the different intrinsic processes involved cannot be stimulated and identified separately. The in silico models are parametric for known components - their behavior is expected to be of a certain shape - and non-parametric for the unknown - as there is no mechanistic reason to give credence to a certain extrapolation. All the in silico components can be simulated separately or in different combinations, allowing the portrayal of different hypothesis. Please click here to view a larger version of this figure.

Figure 2: (A) Coulomb's Law (B) Electric potential field and mobile probe charge (C) Electric current (D) Charge density (E) Current density (F) Equation of continuity (G) Charge conservation law. (A) Electrically charged stationary particles q and Q interact electrostatically through Coulomb's force . (B1) Each charged particle Q generates a scalar field called electric potential at all positions in space: . The maximum work required to move another charged particle q from its position is the product between the charge q and the electric potential generated by Q at position . The electric potential field generated by multiple particles is the sum of the fields generated by each individual particle. (B2) A stationary field with fixed generator particles q and Q, acts with a upon a probe particle with positive charge q_p. In response, q_p moves to minimize its position's electric potential. To describe the motion of q_p, one can derive and the electric field from the electric potential field: . (C) Multiple mobile positively charged probe particles uniformly released in a stationary electric field follow an organized motion. To track the charge configuration without tracking every particle, one can specify at every instant: (D) how space is occupied by particles, assigning a charge density to each infinitesimal volume, according to Gauss's Law, and (E) how particles pass through the boundary surfaces between neighbouring infinitesimal volumes, assigning a current density to each boundary according to Ohm's Law. (F) Charge and current density evolve co-dependently according to the Equation of continuity, as non-uniform particle displacement leads to either accumulation or loss of particles in a certain volume. (G) Within an isolated system, the Charge conservation law prevails and there is no inflow or outflow of charged particles. Notations used:- q,Q,q_p charge and name of the charged particle; - Euclidian norm of the position vector; k - Coulomb's constant; - gradient operator, ε_a - absolute permittivity of medium; σ - conductivity of medium. Please click here to view a larger version of this figure.

Figure 3: Five different levels of complexity for a fibrous mat. SC- smooth with conductive embedded fibers, the simplest model, not accounting for surface morphology or different properties of the constituent components; SNC- smooth with non-conductive embedded fibers; RC- rough with conductive embedded fibers, accounting for surface morphology but not for different component properties; RNC- rough with non-conductive embedded fibers, full proposed model of nanofibrous artificial scaffolds; RNCd- rough with non-conductive embedded fibers coated with two different materials, full proposed model for a sheet of collagen fibers. Length unit used: nanometers. Please click here to view a larger version of this figure.

Figure 4: Example results of the adaptive mesh refinement and the resulting charge density following the simulation. (Left) Automatically generated mesh with extra coarse tetrahedral elements; (Right) Initial mesh adaptively refined during stationary study; smaller elements are required for an accurate result in the areas where simulated signals have abrupt spatial changes. Please click here to view a larger version of this figure.

Figure 5: Fiber alignment angle to electric potential gradient impacts EF strength in surrounding cell culture media when enough complexity is accounted for. SC, SNC, RC, RNC and RNCd are the different levels of complexity for the fibrous mat model presented in Figure 3. Vertical axis marks the alignment angle of the fibers to the electric potential gradient. Abstract electrodes featured - bottom side with high electric potential and top side with low electric potential. Please click here to view a larger version of this figure.

Figure 6: Fiber alignment angle to electric potential gradient impacts space charge density in surrounding cell culture media when enough complexity is accounted for. SC, SNC, RC, RNC and RNCd are the different levels of complexity for the fibrous mat model presented in Figure 3. Vertical axis marks the alignment angle of the fibers to the electric potential gradient. Abstract electrodes featured-bottom side with high electric potential and top side with low electric potential. Please click here to view a larger version of this figure.

Figure 7: Charge movement is influenced by scaffold fiber alignment relative to the EF. Both panels illustrate steady state RNC model predictions. On the left side the fibers are parallel to the EF, while on the right side they are perpendicular. The light red to blue color volume marks charge density, while the arrow volume marks current density orientation. The color of the arrows corresponds to the current density norm. Please click here to view a larger version of this figure.

Name	Expression	Description
Ws	10*Rc*med_ratio	Scaffold width
Ls	10*Rc*med_ratio	Scaffold length
Hs	2*Rf	Scaffold height
med_ratio	1.5	Ratio cell culture media to scaffold
Rc	278.5[nm]	Fiber core radius
r	1.5	Fiber core to coat ratio
Rf	Rc*r	Fiber with coat radius
theta	90[deg]	Fiber orientation angle
Lf	1.3*(Ls*cos(theta)+Ws*sin(theta))	Fiber length
tes	1	Ratio fibre core radius to distance between fibres
n_1	2*(fix((Ws/(2*cos(theta))-Rf)/(2*tes*Rc))+3)*(cos(theta)!=0)+1*(cos(theta)==0)	Max number of fibers if theta<=45
n_2	2*(fix((Ls/(2*sin(theta))-Rf)/(2*tes*Rc))+3)*(sin(theta)!=0)+1*(sin(theta)==0)	Max number of fibers if theta>45
excess	1.2+0.3*abs(sin(2*theta))	First fiber relative offset from scaffold
D	Lf/5	Coat periodicity
prop	0.46	Length of first coat relative to periodicity D
E	100[mV/mm]	Electric field magnitude
V0	E*Ls*med_ratio	Terminal Voltage
omega	500[Hz]	Time dependent study Voltage frequency
p_sigma	0.5	Second coating relative conductivity
p_eps	1.5	Second coating relative dielectric constant

Table 1: Parameters used for simulation

	Culture Media	PEDOT:PSS 1	PEDOT:PSS 2	Collagen Hydrated 1	Collagen Hydrated 2	Silk Fibroin	Collagen Dry
Electrical Conductivity (S/m)	1.7014	1.00E-01	p_sigma * 0.1	2.00E-05	p_sigma * 2e-5	1.00E-08	2.50E-08
Relative Permittivity	80.1	2.2	p_eps * 2.2	9.89	p_eps * 9.89	7.81E+00	4.97

Table 2: Material properties used in simulation

Supplemental Files. Please click here to download this File.

Discussion

The proposed protocol suggests a uniform modelling solution for natural and artificial scaffolds and highlights the need to consider the nanostructure of fibrous scaffolds when inspecting the effects of EF on cells seeded onto such materials. Although a coarse approximation for the EF intensity (electrode potential difference divided by the distance between the electrodes) would lead us to expect a field strength of 100 mV/mm, simulations predict stationary field strengths up to 30% higher in different areas of the mat (Figure 5). This result should be of interest in ES experiment design and data interpretation, as cell death can be caused by too strong EFs. Exposing the electrical microenvironment would enable a direct correlation between ES and cellular development. While several studies present detailed morphology analysis of the used scaffolds^33,43,59, they do not investigate the interplay between the structure, electrical properties of the materials and the EF. This protocol can enable this link, as parameters such as fiber radius, coating layer thickness, distance between fibers and electrical properties of the component materials can be modified according to each experiment by changing the Global Definitions at steps 1.2 and 1.3. Hence, customized 3D spatially resolved charge and current density predictions can be made for both static and dynamic ES regimes.

Scaffold design optimization can be targeted through the RNC and RNCd models with wide parameter range explorations, scaling the proposed morphologies or parts of them. Alternatively, other scaffold configurations can be investigated with the proposed protocol by changing the Array types from Linear to Three-dimensional in section 1.6.5 and adapting Scaffold Geometry in section 1.6.2. However, scaffold optimization cannot be done without an objective. While for tissue engineering purposes the main focus is cell fate, a clearer picture on what stimuli are its main determinants is essential if its reliable control is desired. Charge and current density are good descriptors of cellular electric microenvironments as they show the interplay between the EF and the electrical properties of the different component materials of complex scaffolds such as ECM. The protocol shows how to compute predictions for those metrics given a nanofibrous scaffold geometry and highlights the importance of the alignment angle of the fibers with the EF. Predictions of charge and current density could then be linked to cellular development and thus scaffold and ES regimes may then be optimized for specific tasks.

Interestingly, a study shows that EF exposure generated mechanical stress more than double in strength in composite films with nanofibers perpendicular to the external EF compared to films with parallel alignment⁶⁰. The reported mechanical stress could be a result of Coulomb forces acting between charged fibers, predicted by the rough model simulations (RC, RNC, RNCd) (Figure 6). While these simulations could be useful in investigating this hypothesis, it must be noted that the reported experimental results were obtained in a system with capacitive coupling, and the simulation presents direct coupling.

A limiting factor towards future possible uses of the protocol to estimate a cellular input signal is parameter uncertainty. Geometric uncertain parameters are coating layer thickness and distance between fiber cores. The first one could be inferred by finding the value that leads to a bulk impedance that can be experimentally validated. The second one can be extracted from high resolution material scans. Parameters describing the physical properties of the materials are also affected by uncertainty. However, the electric conductivity and dielectric constant of exemplified materials differ far more than experimental measuring precision (Table 2). Therefore, the reported effects would be maintained despite moderate measurement errors.

The results show how not enough model complexity might hide relevant information. It is important to acknowledge that the protocol simulates a simplified version of the physical phenomenon taking place as it does not account for the different nature of materials involved in the process -conductor (electrodes), semiconductor (coating), dielectric (fiber cores) and electrolytic (surrounding substance) - that are able to influence charge transport. This issue can be accounted for in future model expansions by adding energy transfer delays at the interfaces (i.e., Faradic reactions) and ion transport delays within the electrolyte. Adding complexity should however be guided by experimental validation, as a simple model that reproduces most of what is observed is more useful than a remarkably accurate one that adds little more information but is deeply sensitive to many constituent parameters' uncertainty.

As the end goal of tissue engineering is to create bioreactors that not only mimic one or two aspects of in vivo environments, but replicate and control all cellular developmental cues⁶¹, electromagnetic and mechanical in silico models as well as models of heat transfer between bioreactor components will need to be combined. In a subsequent modelling phase, coupling phenomena between those interactions such as ohmic heating, electrolytic fluid flow, morphological scaffold deformations in response to electrical stimulation⁶⁰ and piezoelectricity⁶² can also be added. However, models should be merged only after each one has been experimentally validated. This way, we can gain a better understanding of each component's influence in the cellular microenvironment, and how stimuli can be optimized.

If the proposed model is experimentally validated, it can be combined with models of biological cells - Figure 1. Charge density patterns and modulations could asymmetrically influence specific ion pumps' activity, impact attachment to the fiber of proteins driving membrane adhesion⁶³ and hence guide migration, proliferation patterns and morphogenesis⁶⁴. Exploring those hypotheses is the way forward in understanding the mechanisms underpinning tissue and cell responses to ES.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Comsol multiphysics 5.2 AC/DC module	COMSOL	-	FEM modelling software