Topological properties of a self-assembled electrical network via ab initio calculation - Scientific Reports

Electrical transportation networks can be found in many disparate areas, including electrical arcs such as lightning^1,2, biological information distribution systems³, the connections between neurons in a brain⁴, and electrical power distribution networks⁵. These type of networks are often not designed or engineered, they grow naturally in accordance to the physical laws that govern their constituents.

Complex flow networks also appear upon careful analysis of other systems. The analysis of complex time series such as EEG data reveals that understanding of the network structure of the generating process is helpful in detecting epileptic seizures⁶. Understanding of the complex network structure of the system dynamics also allows for characterization of oil-water flows^6,7, and gives insight into transitions in nonlinear gas-liquid flows⁸.

Surprisingly, even though the underlying dynamics varies from system to system, certain scaling properties of the resulting networks appear to be universal for a variety of systems⁹. The scaling properties also play an important role in determining the global transportation properties of the network¹⁰. In this work, we consider a system that consists of many electrical conductors which self-assemble into a tree-like network in response to applied electrical voltages or currents¹¹.

Experiments have shown that the degree distribution of the steady state network formed by this system is universal over a wide range of conditions and network sizes. In particular, it was found that the fraction nodes in the network with degree 1 is 25.2% of the total number of nodes, while 23.7% of nodes have degree 3 or higher^11,12. It is interesting to note that this is similar to the degree distribution obtained from other growth processes, notably diffusion limited aggregation¹² and 2D random minimal spanning trees¹³, although it is unknown why this is the case. In addition to the degree distribution, it is also known experimentally¹² that the fractions of nodes with Strahler numbers 1, 2, and 3 is 45.5%, 27.5%, and 16.9% respectively, which is in agreement with Horton’s laws for river networks^14,15. Several global optimization principles have been suggested as governing principles for the self-assembly process, notably a resistance minimizing principle¹⁶ and an entropy production maximization principle¹⁷. The mechanism by which these principles might emerge from the physical dynamics of the self-assembly process still remains unclear.

Some attempts have been made to model the self-assembly process, but these typically involve nonphysical simplifications in order to avoid the complex many-body interactions^18,19,20. These models are unable to predict the scaling properties of the emergent networks, and predict a steady state structure which is qualitatively different from the experimentally observed structure¹⁹. Here we construct a model of the self-assembly process starting from fundamental electrodynamics which includes the many-body interactions by construction. We then develop a method that makes the numerical solution of the model possible. We are then able to calculate the topological properties of the emergent network starting directly from the physical laws of motion. We then use this method to calculate the degree distribution of the network as well as the distribution of nodes with a given Strahler number. This model correctly reproduces the experimentally measured results, and also predicts the topological structure of the emerging network during the formation process. Surprisingly, we find that the observed steady state degree distribution relations are also obeyed during the approach to steady state.

Physical model

We consider the case of i = 1, 2, …, N electrically conducting spheres of radius R. Since the current has been experimentally measured to be quite low²¹, the conductor interactions can be modeled using electrostatics alone. In the empty regions between conductors, the electric potential is determined by the Laplace equation

under the approximation with appropriate boundary conditions.

The total charge on conductor i is related to the potential via

Here the integral is taken over the surface S_i of conductor i and is the surface normal vector of the surface element da_i. The permittivity of the medium is ε.

The conductor charges are taken to be the boundary conditions for Eq. 1 with the exception of an additional conductor i = 0, the ground electrode, which held at a fixed position with zero electric potential. When two or more conductors are in electrical contact, they form a new conductor with a charge equal to the sum of all the conductors in contact. If a conductor is in contact with the ground electrode, it will also have zero electric potential.

The solution of Eq. 1 with these boundary conditions allows the calculation of the electric force on the i^th conductor as another surface integral of the form

This is enough to write down the dynamics of the system:

where m_i, γ_i is the mass and drag coefficient of conductor i, respectively. A charge source term J_s has been added to account for the charge supplied to the system by external means. F(r_i − r_j) is the contact force between conductor i and j. In this work we take the contact force between two spherical objects to be for r_i − r_j < 2 R, and F(r_i − r_j) = 0 otherwise²².

Network analysis

We computed the positions of the N conductors in the system as a function of time for nine values of N between N = 100 and N = 324. These numbers were chosen to be comparable to previous experiments^11,12. A comparison between the steady state produced by the experimental system and the state produced by the numerical solution of the model after t = 120 s can be seen in Fig. 1. This time was chosen such that the system reached an approximately steady state in all cases.

Left: Experimental steady state for N = 512. Right: Numerically calculated state after t = 120 s for N = 289.

Full size image

From the set of conductor positions obtained from the numerical solutions of the equations of motion, Eq. 4, an N node graph was constructed that represents the electrical network at a given time t in which each node corresponds to one conductor. For the analysis, two conductors i, j are considered to be electrically connected at time t if

This is the connection criterion used in the experimental work¹². The weight of the connection between nodes i, j is w_ij = 1 if the conductors are electrically connected and w_i,j = 0 otherwise. We then define the degree d_i of each node i to be

It is also possible to define an anti-arborescence rooted at the ground electrode. To do this, each conductor i is assigned a direct successor D_i with the interpretation that the flow of charge in the network is from i to D_i. The successors are then computed iteratively by defining the successor of all conductors connected to the ground electrode to be the ground electrode. In each subsequent iteration, the successor D_i of the i^th conductor is defined to be the nearest conductor that is connected to i and has a successor provided that i does not already have a successor. This process is iterated until no new successors can be assigned. Depending on the connectivity of , it is possible that not all of the N conductors will be in and so we will use M to denote the number of nodes in .

The total number of nodes which are directly or indirectly connected to the ground electrode and have degree d_i = j, is called Δ_j(t). This is computed as

where δ_i,j is the Kronecker delta.

We define a branch node to be a node i that has a degree d_i ≥ 3. The total number of nodes B(t) which are branch nodes at time t can be calculated as

with the successors D_i defined, each node i can then be assigned a Strahler number s_i. These are determined using the standard procedure. First, we assign s_i = 1 for all with d_i = 1. Then for each node that doesn’t have a Strahler number defined we set s_i equal to the maximum Strahler number of the nodes which have node i as their direct successor. If more than one node with i as its direct successor have the maximum Strahler number, s_i is incremented by s_i → s_i + 1. The number of nodes of with Strahler number j is σ_j:

Comparison to experimental results

Experimental results suggest a linear relationship between and N. Specifically, the relation was found in steady state networks¹¹. Figure 2 shows a plot of Δ₁(t) vs. N computed from the model after t = 120 s of run time compared to experimentally observed steady state relationship.

Bottom: Number of termini Δ₁ as a function of number of conductors connected to the ground electrode M during the formation process out of total N = 324. In both plots the slope of the black line is the experimentally measured value¹¹ 0.252.

Full size image

In addition, the relation was also observed in steady state networks¹¹. Figure 3 shows a plot of B(t) vs. N computed from the model after t = 120 s of run time compared to experimentally observed steady state relationship.

Top: Number of branch points B as a function of total number of conductors N after t = 120 s. Bottom: Number of branch nodes B as a function of number of conductors connected to the ground electrode M during the formation process out of total N = 324. In both plots the slope of the black line the experimentally measured value¹¹ 0.237.

Full size image

A similar result was observed to hold for the distribution of nodes of a given Strahler number¹². Experimental results show that σ₁, σ₂, σ₃ are each linearly related to N. The relations are σ₁ = 0.455 N, σ₂ = 0.275 N and σ₃ = 0.169 N. Figure 4 shows a comparison of the computed distribution of Strahler numbers as compared to the experimentally observed relations.

Black lines are the experimentally observed relations¹².

Full size image

It can be seen from the model that any stationary state of the system must correspond to a connected graph. This is because any conductor that lacks a connection (either directly or through contact with other conductors) to the ground electrode will eventually experience a force directed towards the ground electrode or another conductor in contact with the ground electrode due to the accumulation of charge from the source term J_s. Only in the event that all conductors have a connection to the ground electrode do the forces on the conductors vanish.

Experimentally, it has been noted that the stationary networks rarely have closed loops^11,12. This may be due to the form of Eq. 3, which shows that the force on any surface element of a conductor is normal to the surface and directed outwards. A closed loop can be thought of a single conductor with electric field inside the loop. Any such loop may experience a force that acts to expand the loop, and thus separate the conductors that comprise it. This force can only be zero in the case that everywhere on the outside surface of the loop as well. Therefore, closed loops in the conductor network are at best neutrally stable, and unstable in the presence of any external electric field.

The stationary states of the system are thus expected to be connected acyclic graphs, or trees. For any M node tree , we can define the fraction f_i(t) of nodes in that have degree i as

The fractions f_i are not completely independent. First there is a normalization constraint with maximum degree i_max

Second, there is a constraint related to the number of branches and endpoints. This can be interpreted as a statement that every branch in the tree ultimately has an end associated with it, and the number of additional branches created by a node of degree i is i − 2. This can be expressed in terms of the fractions f_i as

In the limit , the 1/M term may be dropped, and equations 11 and 12 combine to give

Eq. 13 implies that the average degree of nodes is 〈i〉 = 2 for large M.

For spherical particles in 2D, the maximum degree any node can have i_max = 4, as any degree larger than this would be a closed loop if connections are defined using 5. Thus there are only two independent numbers that specify the degree distribution of the network. The degree distribution can be expressed in terms of the fraction of endpoints f₁ and the fraction of branch nodes b = f₃ + f₄ as follows.

Using the experimental values of f₁ = 0.252 and b = 0.237 we obtain the full degree distribution

We note the apparent similarity of this degree distribution with that minimum spanning trees for random sets of nodes¹³ which have f₁ = 0.253, f₂ = 0.527, f₃ = 0.204, and f₄ = 0.016. This may suggest that the network evolves in accordance with some optimization principle, such as the resistance minimizing principle¹⁶, or the maximum entropy production principle¹⁷.

The model is able to reproduce the experimental degree distribution for self assembled electrical networks starting from the physical interaction between conductors as well as correctly explain the rarity of closed loop structures. In addition, the model predicts that the degree distribution of the network remains constant during the approach to the steady state. We are also able to reproduce the experimentally observed distribution of Strahler numbers in the network.

Calculation of electric potential

Solutions to Eq. 1 were generated in parallel with the red-black Gauss Seidel method²³. The region of interest was a 2D square L = 5.12 cm on each side which was discretized into a square grid with 512 × 512 sites. Neumann boundary conditions were taken at the edges of the square region, and all sites within conductor i = 0 were held at zero potential. In addition, a constraint was imposed such that the total charge on each conductor with i > 0 remained fixed at it’s known value from Eq. 4. This constraint corresponds to finding a set of potentials {ϕ_i} on the remaining conducting surfaces which were used as Dirichlet boundary conditions.

This is accomplished by first defining the function

In this equation is the solution of Eq. 1 subject to the Dirichlet boundary conditions given by the set of ϕ_i’s. Equation 16 gives the charge on the i^th conductor if all the conductors were held at the known fixed voltages {ϕ_i}. With this, we construct the functions

Along with the total error associated with these initial conditions

By the uniqueness property of solutions to Eq. 1, there is only one solution of the form E({ϕ_i}) = 0 for Eq. 18. Therefore this solution gives the correct conductor potentials {ϕ_i} corresponding to the charges {Q_i}.

The problem of solving for the conductor voltages is now a root finding problem which can be solved by gradient descent. Gradient descent for this problem would involve making the update

for some positive γ_i. However, this update requires N sums over N terms, and also requires knowledge of the coefficients c_ij. Instead, we make a local update by considering only one term of the sum

Here we have made the redefinition γ_ic_ii → γ since c_ii is always positive as well^24,25. This update does not require the sum over N terms, and does not explicitly require knowledge of c_ii. It only requires that γ be chosen to ensure convergence. At each iteration, the conductor potentials are changed by an amount

From now on we will write . This update causes the errors e_i to change by

In general this can be either positive or negative. The total error E then changes on the (n + 1)^th iteration by

Since c_ij is a positive matrix^24,25, the quantity is always positive. Therefore, the total error is always decreasing provided γ is chosen such that

Thus the update rule in Eq. 20 converges provided γ is chosen to satisfy Eq. 24. In this work a value of γ = 7.95 nV/C was found to be sufficient for convergence. This method was iterated until a solution with E/N < 2.3 pC was obtained.

Integration of the equations of motion

From the solution to Eq. 1 on the discretization grid, a bilinear interpolation was used to compute the potential between grid points. This is given by

for x, . Here x is in units of grid spacing (px). This interpolation allows the computation of the electric potential in the area between grid points by approximating the behavior in the local region as changing linearly in x and y between the known values on the grid points. The surface integral in Eq. 3 was calculated from this interpolation using the trapezoid method of integration with N_s = 500 evenly spaced points around each conductor.

where x_ij = x_i + R cos (2πj/N_s), and y_ij = y_i + R sin (2πj/N_s). Here the total force vector is approximated as the vector sum of the electrical force exerted on each of the N_s discrete surface elements. Similarly the conductor charges from Eq. 2 are computed from

Here the charge is computed as the scalar sum of the charge on each of the N_s discrete surface elements.

The equations of motion in Eq. 4 were numerically integrated using Euler’s method with a timestep of Δt_m = 0.01 ms for the mechanical degrees of freedom and Δt_e = 0.001 s for the electrical degrees of freedom.

The mass of all conductors was set to m_i = 16 g. The fluid drag was assumed to be Stokes drag, and so γ_i = 6πμR with μ = 650 cP for castor oil²⁶. The permittivity of the medium was set to be ε = 4.7ε₀, where ε₀ is the permittivity of free space. This is near the dielectric strength of castor oil²⁷. For the contact dynamics we used . Conductors were considered to be in contact if |r_i − r_j| <= 2 R + 2 px. The charge source was J_s = 3.8 pC/s.

The conductors were initially laid out in a square grid centered in the region of interest with a separation between their centers of 3 R and the radius of all conductors was set to be R = 0.8 mm. The ground electrode position was fixed at (x₀, y₀) = (496, 257). The position of the ground electrode is shifted one grid site off of the center line in the y direction to explicitly break symmetry.

How to cite this article: Stephenson, C. et al. Topological properties of a self-assembled electrical network viaab initio calculation. Sci. Rep. 7, 41621; doi: 10.1038/srep41621 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

1. Yoshida, S., Akita, M., Morimoto, T., Ushio, T. & Kawasaki, Z. Propagation characteristics of lightning stepped leaders developing in charge regions and descending out of charge regions. Atmospheric Research 106, 86–92 (2012).

   Article  ADS  Google Scholar 

2. Shi, W., Li, Q. & Zhang, L. A stepped leader model for lightning including charge distribution in branched channels. Journal of Applied Physics 116, (2014).

3. Cai, L., Tabata, H. & Kawai, T. Self-assembled dna networks and their electrical conductivity. Applied Physics Letters 77, (2000).

4. Buzsàki, G. & Draguhn, A. Neuronal oscillations in cortical networks. Science 304, 1926–1929 (2004).

   Article  ADS  Google Scholar 

5. Rohden, M., Sorge, A., Witthaut, D. & Timme, M. Impact of network topology on synchrony of oscillatory power grids. Chaos 24 (2014).

6. Gao, Z. K., Cai, Q., Yang, Y. X., Dang, W. D. & Zhang, S. S. Multiscale limited penetrable horizontal visibility graph for analyzing nonlinear time series. Scientific Reports 6, 35622 (2016).

   Article  CAS  ADS  Google Scholar 

7. Gao, Z. K., Yang, Y. X., Fang, P. C., Jin, N. D., Xia, C. Y. & Hu, L. D. Multi-frequency complex network from time series for uncovering oil-water flow structure. Scientific Reports 5, 8222 (2015).

   Article  CAS  Google Scholar 

8. Gao, Z. K., Fang, P. C., Ding, M. S. & Jin, N. D. Multivariate weighted complex network analysis for characterizing nonlinear dynamic behavior in two-phase flow. Experimental Thermal and Fluid Science 60, 157–164 (2015).

   Article  Google Scholar 

9. Barabàsi, A. L., Dezso, Z., Ravasz, E., Yook, S. H. & Oltvai, Z. Scale-free and hierarchical structures in complex networks. AIP Conference Proceedings 661 (2003).

10. Xu, P. & Yu, B. The scaling laws of transport properties for fractal-like tree networks. Journal of Applied Physics 100 (2006).

11. Jun, J. K. & Hübler, A. H. Formation and structure of ramified charge transportation networks in an electromechanical system. Proceedings of the National Academy of Sciences of the United States of America 102, 536–540 (2005).

    Article  CAS  ADS  Google Scholar 

12. Soni, V. H., Ketisch, P. M., Rodrìguez, J. D., Shpunt, A. & Hübler, A. W. Topological similarities in electrical and hydrological drainage networks. Journal of Applied Physics 109, 036103–036103 (2011).

    Article  ADS  Google Scholar 

13. Ketisch, P. M., Rodrìguez, J. D. & Hübler, A. W. Modeling the degree distribution of a fractal transportation network with a minimum spanning tree graph. Santa Fe Institute http://santafe.edu/media/cms_page_media/263/JuanDiegoRodriguez.pdf (2008).

14. Horton, R. E. Erosional development of streams and their drainage basins; hydrophysical approach to quantitative morphology. Geological society of America bulletin 56, 275–370 (1945).

    Article  ADS  Google Scholar 

15. Dodds, P. S. & Rothman, D. H. Unified view of scaling laws for river networks. Phys. Rev. E 59, 4865–4877 (1999).

    Article  CAS  ADS  Google Scholar 

16. Dueweke, M., Dierker, U. & Hübler, A. Self-assembling electrical connections based on the principle of minimum resistance. Phys. Rev. E 54, 496–506 (1996).

    Article  CAS  ADS  Google Scholar 

17. Belkin, A., Hubler, A. & Bezryadin, A. Self-assembled wiggling nano-structures and the principle of maximum entropy production. Scientific Reports 5, 832 (2015).

    Article  Google Scholar 

18. Hadwich, G., Mertè, B., Hübler, A. & Lüscher, E. Stationary dendritic structures in an electric field. Helvetica Physica Acta 63, 487–488 (1990).

    Google Scholar 

19. Marani, M., Banavar, J. R., Caldarelli, G., Maritan, A. & Rinaldo, A. Stationary self-organized fractal structures in an open, dissipative electrical system. Journal of Physics A: Mathematical and General 31, L337 (1998).

    Article  ADS  Google Scholar 

20. Mertè’, B. et al. Formation of self-similar dendritic patterns with extremal properties. Helv. Phys. Acta 62, 294–97 (1989).

    Google Scholar 

21. Stephenson, C. & Hubler, A. Stability and conductivity of self assembled wires in a transverse electric field. Scientific reports 5, 15044 (2015).

    Article  CAS  ADS  Google Scholar 

22. Hanaor, D. A., Gan, Y. & Einav, I. Contact mechanics of fractal surfaces by spline assisted discretisation. International Journal of Solids and Structures 59, 121–131 (2015).

    Article  Google Scholar 

23. Gourlay, A. R. & McGuire, G. R. General hopscotch algorithm for the numerical solution of partial differential equations. IMA Journal of Applied Mathematics 7, 216–227 (1971).

    Article  MathSciNet  Google Scholar 

24. Diaz, R. A. & Herrera, W. J. The positivity and other properties of the matrix of capacitance: physical and mathematical implications. Journal of Electrostatics 69, 587–595 (2011).

    Article  Google Scholar 

25. Herrera, W. J. & Diaz, R. A. The geometrical nature and some properties of the capacitance coefficients based on Laplace’s equation. American Journal of Physics 76, 55–59 (2008).

    Article  ADS  Google Scholar 

26. Quinchia, L., Delgado, M., Valencia, C., Franco, J. & Gallegos, C. Viscosity modification of different vegetable oils with EVA copolymer for lubricant applications. Industrial Crops and Products 32, 607–612 (2012).

    Article  Google Scholar 

27. Mutlu, H. & Meier, M. A. R. Castor oil as a renewable resource for the chemical industry. European Journal of Lipid Science and Technology 112, 10–30 (2010).

    Article  CAS  Google Scholar 

Download references

Authors and Affiliations

1. Department of Physics, University of Illinois at Urbana Champaign, Urbana, Illinois, USA

   C. Stephenson, D. Lyon & A. Hübler

Contributions

C.S. is the lead author and performed the numerical simulations and analysis. C.S., D.L., and A.H. conceived the model. A.H. supervised the work. All authors generated ideas for this paper and discussed the text at all stages.

Corresponding author

Correspondence to A. Hübler.

Competing interests

The authors declare no competing financial interests.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Reprints and permissions