Applications to Biological Networks of Adaptive Hagen-Poiseuille Flow on Graphs

Ana Filipa Martinho Valente

Engineering Physics

:     Rui Manuel Agostinho Dilão

:        Luís Humberto Viseu Melo

:         Rui Manuel Agostinho Dilão

Member of the Committee:	    Ana Maria Ribeiro Ferreira Nunes

June 2023

"No matter where you are,

everyone is always connected."

I would first like to express my heartfelt gratitude to my supervisor, Professor Rui Dilão, for his insurmountable patience and availability, for the compelling and thought-provoking conversations, for the rigorous counseling and invaluable feedback, for the kindness, compassion and encouragement he has shown me, and for all the expertise and wisdom he so kindly shared with me. Thanks to his continuous support, to everything he has taught me, and to his commitment to quality academic work, I have a newfound love for research. His guidance has most definitely not only enriched my academic journey but also encouraged my personal growth.

I'm forever indebted to my wonderful dear boyfriend, André Pereira, for all the unwavering love and support he has shown me throughout not only this difficult, jarring time but ever since we first met. Thank you so much for all your gentle and kindhearted words, for the endless motivation you provide, for all the time we spent relaxing, going on adventures and celebrating, and for helping me flourish into a person I can be ever more proud of. Being able to love you is an enormous privilege.

I would like to deeply thank all my dearest friends, in particular Inês Ferreira, Rita Silva, Bernardo Barbosa, Tomás Lopes and José Maria Nina Carreira Ferreira da Cruz. The conversations we had, the games we played, the music we shared, the cats and the funny images we laughed at together during this intense time all helped me immensely to stay sane and energetic enough to fight through life. I wish you all the very best, peaceful, happy times in your life.

Last but not least, I'm eternally grateful for my family, in particular my parents, for all their unconditional support and caring questions. Despite having to go through some of the toughest times in their life, they have shown me nothing but love and patience and support throughout these last few years, even when I proved difficult and stubborn.

I appreciate and cherish every single one of you. I wish you all the very best in your life.

I'd like to extend my gratitude to Georgios Cherouvim for allowing his mesmerizing artwork to be displayed on the cover page of this thesis. Cover page image adapted from <cit.>.

CHAPTER: RESUMO

Physarum polycephalum é um protista unicelular com numerosos núcleos cujo corpo é constituído por uma rede de veias. À medida que explora o seu meio, adapta-se e otimiza a sua rede, tendo em conta estímulos externos. Foi demonstrado que exibe comportamentos complexos, como resolver labirintos, encontrar o caminho mais curto e criar redes robustas, eficientes e económicas. Vários modelos foram desenvolvidos para tentar simular a adaptação da sua rede para compreender os mecanismos por detrás do seu comportamento e desenvolver redes eficientes.
Esta tese pretende estudar um modelo recentemente desenvolvido e fisicamente consistente baseado em fluxos de Hagen-Poiseuille adaptativos; para tal, irá determinar propriedades das árvores produzidas pelo modelo e irá examiná-las para determinar se são realistas e consistentes com a experiência. Esta tese também pretende usar o mesmo modelo para produzir redes curtas e eficientes, aplicando-o a uma rede de transporte real.
Observámos que o modelo é capaz de criar redes que são consistentes com outras redes biológicas: seguem a lei de Murray no estado estacionário, mostram estruturas semelhantes às presentes nas redes do Physarum e ainda exibem peristalse (oscilação dos raios das veias) e shuttle streaming (o movimento de trás para a frente do citoplasma do Physarum) em algumas partes das redes.
Usámos também o modelo em conjunto com diferentes algoritmos estocásticos para produzir redes curtas e eficientes; quando comparadas com a rede ferroviária de Portugal continental, todos os algoritmos produziram redes mais eficientes que a rede real e alguns produziram redes com melhor relação custo-benefício.

Palavras-chave:     Physarum polycephalum; Fluxo de Hagen-Poiseuille; Rede adaptativa;

Optimisação de redes; Eficiência de transporte; Árvore mínima de Steiner

CHAPTER: ABSTRACT

Physarum polycephalum is a single-celled, multi-nucleated slime mold whose body constitutes a network of veins. As it explores its environment, it adapts and optimizes its network to external stimuli. It has been shown to exhibit complex behavior, like solving mazes, finding the shortest path, and creating cost-efficient and robust networks.
Several models have been developed to attempt to mimic its network's adaptation in order to try to understand the mechanisms behind its behavior as well as to be able to create efficient networks.
This thesis aims to study a recently developed, physically-consistent model based on adaptive Hagen-Poiseuille flows on graphs, determining the properties of the trees it creates and probing them to understand if they are realistic and consistent with experiment.
It also intends to use said model to produce short and efficient networks, applying it to a real-life transport network example.
We have found that the model is able to create networks that are consistent with biological networks: they follow Murray's law at steady state, exhibit structures similar to Physarum's networks, and even present peristalsis (oscillations of the vein radii) and shuttle streaming (the back-and-forth movement of cytoplasm inside Physarum's veins) in some parts of the networks.
We have also used the model paired with different stochastic algorithms to produce efficient, short, and cost-efficient networks; when compared to a real transport network, mainland Portugal's railway system, all algorithms proved to be more efficient and some proved to be more cost-efficient.

Keywords: 	    Physarum polycephalum; Hagen-Poiseuille flow; Adaptive network;

Network optimisation; Transport efficiency; Steiner minimal tree

[r]$(i,j)$Edge connecting nodes $i$ and $j$.

[r]$L_{ij}$Length of edge $(i,j)$.

[r]$L$Total length of graph tree.

[r]$r_{ij}$Radius of edge $(i,j)$.

[r]$Q_{ij}$Volumetric flux flowing through edge $(i,j)$.

[r]$D_{ij}$Conductivity of edge $(i,j)$.

[r]$p_{i}$Fluid pressure at node $i$.

[r]$p_{ij}$Pressure difference between nodes $i$ and $j$.

[r]$\mathcal{G}$Graph (set of nodes connected by a set of edges).

[r]$\mathcal{V}$Set of nodes of a graph.

[r]$E$Set of edges of a graph.

[r]$x_i^n$Spatial coordinate of node $i$ along the $n$th spatial direction.

[r]$\eta$Dynamic viscosity of fluid.

[r]$S_i$Inward/outward volumetric flux of fluid to/from the network at node $i$, making node $i$ a source/sink, if $S_i > 0$/$S_i < 0$, respectively.

[r]$V_{ij}$Volume of edge $(i,j)$.

[r]$V$Total volume of fluid of the network.

[r]$\mathcal{P}$Total dissipated energy per unit time of the network.

[r]$\tau$Re-scaled time variable: $ct$, with $c$ a positive constant.

[r]$I_0$Total inward volumetric flux given by all the source nodes to the network.

[r]$\nu$Average length between all possible pairs of sites.

[r]$E'$Set of effectively conducting edges of a graph (edges whose conductivity value is above a certain threshold).

[r]$L_{s_i s_j}$Length of shortest path in tree between sites $s_i$ and $s_j$.

[r]$\sigma$Standard deviation.

[r]$\overline{a}$Mean value of variable $a$.

[r]$\Delta\tau$Time step used in the adaptive conductivities network algorithm.

[r]$P$Perimeter of a polygon.

[r]$L_{\text{Stei}}$Theoretical length for the minimum Steiner tree for a polygon.

[s]$i,j,k,m$Indices pertaining to nodes (single index) or edges (double index) of the network.

[type=, title=]

PhysarumPhysarumPhysarum polycephalum.
FSFSFood source(s).
SMTSMTSteiner minimum tree.
MSTMSTMinimum spanning tree.
FTFTFault tolerance.
TLTLTotal length.
Freq. inv.Freq. inv.Frequency of inversion.
Rel. str.Rel. str.Relative strength.

CHAPTER: INTRODUCTION

§ MOTIVATION

There is a multitude of dynamical systems (biological, physical, electrical) that can be described using networks and graphs. Among them, transport networks are prevalent. Transport networks are essential, as they must allocate information and resources as efficiently and quickly as possible throughout a system. Some examples of transport networks whose study remains relevant today are the phloem and xylem of plants, the network of channels that transport nutrients in their leaves <cit.>, the vein irrigation network of a tumor <cit.>, human transit networks, such as road or railway systems <cit.>, various electrical networks <cit.> and organisms like Physarum polycephalum, more commonly referred to as slime mold (described in section <ref>), or Dictyostelium discoideum <cit.>, which display highly adaptive vessel networks. There are also a number of mathematical problems and formal systems that constitute relevant and difficult network optimisation problems, such as Steiner tree geometries <cit.>, the traveling salesman problem <cit.> and the first passage percolation <cit.>, whose results still have important applications in physics, engineering and biology today. The study of these real-world networks and these theoretical problems can lead to results that are often applicable to many other different transport networks and systems.

Key aspects of these networks can be described using hydrodynamics. While some of the transport networks do in fact transport a fluid of some sort, such as the cardiovascular system, the xylem and phloem, or the adaptive network of veins of organisms like Physarum polycephalum, many other transport networks can be approximately described using fluid dynamics. Specifically, many relevant systems can be described by the Hagen-Poiseuille formalism (which is described in section <ref>), which describes the laminar flow of viscous fluids through a cylindrical vessel of constant cross section. Some biological systems of relevance that can be described this way are blood vascular systems and Physarum polycephalum.

Physarum polycephalum specifically has recently become a prolific topic of discussion, as this acellular protist displays high-level behaviours, such as solving mazes <cit.>, despite not having centralized control or a neurological system. The network formed by its veins adapts dynamically to its environment; this adaptation consists of the formation of new veins, the destruction of old, less important paths, and the modification of vein radii over time. Due to these adaptations, in the presence of several different food sources it can create networks with a comparable efficiency, fault tolerance and cost to real human-made networks (such as the Tokyo rail network) <cit.>. By studying this brainless, simple organism, one could understand the mechanisms behind its obtainment and processing of information and decision-making. After all, Physarum and all living beings have gone through the process of evolution for millions of years, and have thus optimised their transport networks greatly; they can provide great inspiration for the development of efficient problem-solving algorithms.

Modeling this type of Hagen-Poiseuille adaptive flow in these biological networks is important and can provide knowledge about many diverse important topics. Namely, it may grant more understanding regarding efficient, robust and cost-effective network formation, growth and adaptation. This understanding could then be applied to several different areas and challenging problems; for example, it could provide insight into angiogenic processes like cancer development.

§ OBJECTIVES

This thesis has two main intentions: first, to study the adaptive Hagen-Poiseuille flows model, analyzing several properties and consequences of the adaptive flow on graph structures and flow behavior, and understand how reliable it is at predicting real physical results; second, to apply said model to a specific graph theory problem and determine how effective the model is at solving it.

The first objective comprises several additional goals. Said goals are 1) determining how the model behaves when the mesh it is applied to reaches the continuous limit, 2) establishing whether or not Murray's law is verified dynamically for this model, 3) studying properties of simple networks and how they compare to real results, 4) comparing the results of the model to the patterns exhibited by Physarum polycephalum and 5) ascertaining whether or not the model exhibits key physical phenomena shown by Physarum: shuttle streaming and peristalsis.

The second objective, which essentially is utilizing the model to find contrasting trees of different efficiencies and extremal lengths, will be achieved by applying the model to two different types of configurations (geometric regular polygons and a real-world transport network case) and using dynamic algorithms to obtain said different trees.

§ THESIS OUTLINE

This thesis is organized into five chapters. Chapter <ref> starts by introducing the work's subject matter and objectives and the motivation behind the thesis.

Chapter <ref> describes into detail Physarum polycephalum's physical characteristics and its relevant intelligent behaviors, as well as the state-of-the-art network models used to mimic its behavior. It also introduces and explains the adaptive Hagen-Poiseuille flows model (that aims to describe Physarum's network dynamics) that will be the basis of the work of this thesis, and this chapter concludes with some notes regarding some graph theory problems that will be addressed later.

Chapter <ref> delineates the algorithms and methods used to put said adaptive flows model into practice and tests the model regarding some relevant physical properties that are observed in the simulations. It concludes with a direct application of the model to the comparison of obtained trees with Physarum's real network attributes.

Chapter <ref> illustrates another application of the model in question: providing stochastic solutions to the Steiner minimal tree problem in graph theory. The stochastic algorithms used to create the steady-state trees and the parameters used to evaluate them are first introduced and are later applied to two different cases: first, to simple geometric polygons, and then, to a real communication network (the case of mainland Portugal).

To complete the work, chapter <ref> provides a broad synoptic analysis of the results presented in this thesis and the conclusions that are possible to derive from them, as well as suggests several additional work topics and improvements to the model that could derive from this work.

CHAPTER: PHYSARUM POLYCEPHALUM AND ITS NETWORK DYNAMICS

§ PHYSARUM POLYCEPHALUM

Physarum polycephalum, also known as slime mold or "the blob", has been the subject of many biological and biochemical studies throughout the 20th century and, in recent years, of physical studies. Some of the main goals of the modern studies are to shine a light on mechanisms of information processing, along with studying self-regulated networks and their topology.

Physarum has been studied for many reasons, namely for its simplicity regarding culturing and handling in the lab, its macroscopic size (which facilitates observation), and its unique "intelligent" behavior.

I will start by describing slime mold physically in section <ref>; then, I will discuss certain relevant oscillatory phenomena theorized to be at the source of Physarum's motion and intelligent behavior in section <ref>; finally, I will describe said intelligent behavior that Physarum is reported of showing in section <ref> and discuss its implications.

§.§ Physical characteristics

Physarum polycephalum is a unicellular amoeboid protist with multiple nuclei.

Despite being called "slime mold", it is not a fungus, but rather it displays characteristics shared by fungi, animals and plants. Its complex life cycle allows it to have strong versatility and survivability; for example, Physarum can enter a dormant state for up to several years in non-favorable environmental conditions, and said state can be reversed easily by providing it with water and nutrients <cit.>.

Physarum polycephalum's plasmodium in its natural environment, growing on tree bark. a) An amorphous yellow mass forms after some time of growth. b) The food foraging and growth process for Physarum involves the creation of a network of veins of protoplasm, which adapts over time. Images obtained from a) <cit.> and b) <cit.>.

The most frequently found form of Physarum is the plasmodium. In the plasmodial stage of its life cycle, it is a single extremely large multinucleate cell containing millions of nuclei. It appears as a macroscopic, bright yellow amorphous mass that can grow up to a size of 10 m$^2$. It prefers to live in damp and dim habitats and is often seen growing on the sides of trees (see figure <ref>) <cit.>.

It feeds on various organic materials, such as bacteria, mushrooms, fungal spores, and decaying matter (as well as cornflakes when fed in a laboratory environment), by covering the food with its protoplasm, which is later digested in the body by enzymes. To procure food sources, this living being can move at speeds of about 1-4 cm/h <cit.>; it spreads its protoplasm radially outward, in the shape of a network of tubular veins. As Physarum grows, its network's topology is modified and optimised (see figure <ref>).

Physarum polycephalum's plasmodium in a laboratory environment. a) Notice the network created between different food sources (the yellow blobs): a single, emphasized vein forms between the food sources, surrounded by a number of tiny, thin veins; at the top of the image, Physarum's network grows indiscriminately in the direction it has yet to explore. b) Physarum's new growth is signaled by numerous, thin veins growing outwards in fan-like patterns (see the edges of the network), while the older part of the network shows optimised, few, thick veins connecting to the edges of the network. Images obtained from a) <cit.> and b) <cit.>.

Physarum's cell is not only made up of an extremely large number of nuclei, but also of different forms of cytoplasm (that is, the fluid inside the cell that is located outside of the nuclei): endoplasm, which is the fluid cytoplasm that moves throughout the cell, and ectoplasm, which is the gel-like rigid cytoplasm that makes up the outer membrane of the cell and contains the cytoskeleton of the cell. This cytoskeleton is made up of a system of proteins, notably actin and myosin. While the filaments of actin are responsible for the structural support of the vein walls, myosin is the motor protein involved. <cit.>

§.§ Dynamic, periodic phenomena

A solid grasp on Physarum's locomotion is necessary to construct a proper toy model of this organism. Physarum polycephalum locomotion is ensured by two periodic phenomena: shuttle streaming and peristalsis, and these two are theorized to be connected.

Peristalsis is a wave of cross-sectional contractions across a tubular vein of the network. The interactions between the proteins that make up the cytoskeleton of the cell, actin and myosin, create relaxation-contraction cycles of the vein walls, which propels the movement of the endoplasm.

Shuttle streaming describes the movement of the endoplasm through the network veins. Due to peristalsis, the fluid moves in a back-and-forth pattern throughout the veins of the entire organism (with a period of about 100s). This allows for the distribution of resources throughout an entire cell.

These rhythmic oscillations and their timely coordination are what allow the organism to move in a certain direction. A gradient of pressure is produced by these waves, which propels the movement of endoplasm towards the leading edges (where the growth of ectoplasmic proteins occurs at the same time). <cit.>

Physarum individuals are also able to maximize their internal endoplasmic ﬂows by adapting these contraction waves to their size, thus optimizing transport.
It was shown that the transport is optimal when the wavelength of the peristaltic wave is of the order of magnitude of the size of the network. Thus, Physarum, despite not having a nervous system, is able to coordinate its growth and adapt the flow to its size.
These adaptive patterns come about due to the interaction of a global mass conservation constraint with a local constraint: minimizing the phase difference between neighboring vein sections.
These contractions are also driven by external stimuli, as the amplitude and frequency of the contractions are altered across the whole organism. Attractive stimuli, like food sources, usually increase the amplitude and frequency of the contractions, while repulsive stimuli, like excessive lighting, decrease said amplitude and frequency. <cit.>
Even though it is not yet well understood how all these oscillations are coordinated throughout the entire organism and how they arise from stimuli, it is theorized that a signaling molecule could be involved: this molecule could cause local increases in contraction amplitude as it was transported along with the endoplasm and could initiate a feedback loop to propel its own motion <cit.>.

These biochemical and hydrodynamic oscillations thus allow Physarum to grow its network and to dynamically adapt its morphology to external stimuli.

§.§ Intelligent behaviour observed

Even though Physarum is a comparatively simple protistic organism, with no nervous system or centralized control, it has been shown to exhibit high-level intelligent behavior and create highly optimized networks. These behaviors are now described.

A famous experiment was conducted to determine whether or not Physarum was able to solve a maze, and the protist did so successfully <cit.>. A Physarum specimen was made to cover an entire maze with agar substrate and plastic film walls (agar acts as food for the specimen, while the plastic film is a dry surface with the specimen avoids). After an operator placed food (agar blocks) at the entrance and exit of the maze, Physarum rearranged itself: it eliminated useless veins and enlarged the most favorable ones. It created the most efficient and shortest path between the two food sources every time, thus solving the maze (figure <ref>).

Later experiments had Physarum recreate different man-made transportation systems; most notably, Physarum recreated the Tokyo region's railway system. Food sources (oat flakes) were placed in such a way that they mimicked the location of major cities in the region and a Physarum specimen was placed on a central food source. Topographical limitations like mountainous terrain or bodies of water were also taken into account, and were replicated using light sources, as Physarum avoids luminous locations. As Physarum grew out of the central location (radially outward) and it found the food sites surrounding it, it adjusted its veins and the morphology of its network (figure <ref>). This network was found to have "comparable efficiency, fault tolerance and cost" to the Tokyo railroad network. Note that "cost" refers to the total length of the network, "efficiency" means the average minimum distance between each pair of food sources, and "fault tolerance" means the probability of disconnecting a part of the network by removing a certain link. <cit.>

These experiments showcase one of the important applications of studying Physarum. Designing transport networks that take into account both robustness and cost-efficiency is a modern, ever-present problem. Physarum is able to create such networks naturally as part of its growth process, and efficient algorithms can be developed inspired by it.

Physarum polycephalum's intelligent behaviour experiments. a) Physarum solving a maze over time. First, the specimen is made to cover the entire maze and food sources are placed at each end of the maze (left); over time, it adapts its network (middle) until it reaches the shortest length network possible connecting both food sources/ends of the maze, which is the solution to the maze (right). b) Physarum constructing a network resembling the Tokyo rail system over time. The specimen is placed in a central location, while food sources are placed in the geographical location of major cities of the Tokyo region (A - 0hr); Physarum starts by growing radially from its original location (C - 8hr) until it reaches all food sources (E - 16hr); as it grows, it adapts its network, eliminating unnecessary links and reinforcing central ones, creating a cost-efficient and robust network (F - 26hr). Images adapted from a) <cit.> and b) <cit.>.

There are even more intelligent behaviours associated with Physarum. For instance, slime mold has also been found to exhibit memory, both internal <cit.> and external <cit.>. As Physarum moves around its environment, it leaves behind a trail of mucus and, when the organism comes across its trail on its search for sustenance, it avoids the already explored areas - this is what one means by external memory. Additionally, Physarum can also correctly choose the most appropriate and efficient diet for itself when provided with several options <cit.>.

Even though several protistic or fungic species (in addition to Physarum polycephalum) show intelligent behaviour despite not having a nervous system, the mechanisms behind this intelligence are not yet comprehended. The goal of studying this brainless, simple organism is to understand the mechanisms behind the obtainment and processing of information and decision-making.

§ PHYSARUM NETWORK MODELS

Physarum polycephalum's plasmodium stage is quite complex, and there exists no one model yet to describe it fully. Models that describe Physarum usually focus on one of three connected aspects that are ultimately responsible for Physarum's behaviour: its network's adaptation, its growth or the oscillatory patterns described in section <ref>. This thesis shall focus on Physarum's network adaptation, and on a particular recently-developed model.

Physarum's network optimization models are used to not only better understand the mechanisms behind this organism's behavior, but to also solve different graph theory problems or related real-world issues <cit.>.

There are two main models described in the literature: the multi-agent model proposed by Jones <cit.> and the Physarum Solver proposed by Tero et al. <cit.>.

The multi-agent system is a phenomenological model that states that the macroscopic complex network modifications observed in Physarum networks are a product of simple microscopic interactions between small parts of the plasmodium that acted as particle-like agents. While this model was able to obtain shortest path-like patterns and Steiner minimum tree-like solutions, it is a biologically unrealistic model as agents do not act like a continuous network and always move towards areas of positive stimuli, which is not accurate to Physarum's food searching process.

The Physarum solver, on the other hand, is a model based on the hydrodynamics of the endoplasm, in which network optimization happens through a feedback loop between the flows of the network and the size of the tubular veins that make it up. This model proved to be more biologically realistic than the multi-agent system and it showed resulting networks similar to those created by Physarum. However, the model isn't physically consistent, as it does not conserve the total volume of fluid of the network.

The model used to describe Physarum's network optimization dynamics in this thesis was developed by Almeida and Dilão <cit.>. It is based on Tero's Physarum Solver but solves its biggest physical inconsistency: the fact that volume needs to be conserved. It also uses dissipated energy minimization considerations to obtain its adaptation equations. It is described in the following section <ref>.

§ ADAPTIVE HAGEN-POISEUILLE FLOW ON GRAPHS

In this section, the model derived by Almeida and Dilão <cit.> used in this thesis to illustrate the Physarum network optimization process is described. First, an introduction to the Hagen-Poiseuille law, which is essentially the physical basis of the model, is done in section <ref>. Then, an overview of the derivation of the model is done in section <ref>. For a full derivation, please refer to <cit.>. Finally, a biologically relevant law illustrated in this model called Murray's law is discussed in section <ref>.

§.§ Hagen-Poiseuille flow

The Hagen-Poiseuille (H-P) flow refers to the steady laminar flow of an incompressible, viscous fluid (that is, a fluid with a low Reynolds number) through a cylindrical channel with slippery boundary conditions.
The formalism is obtained by applying several simplifications to the Navier-Stokes equation for a cylindrical channel, namely symmetry considerations, assuming steady-state conditions and incompressible fluid considerations, and also taking into account mass conservation.
It is often used to describe blood flow <cit.> and Physarum polycephalum's endoplasm flow <cit.>.
For a channel $(i,j)$ (that is, between points/nodes $i$ and $j$), with length $L_{ij}$ and radius $r_{ij}$, the Hagen-Poiseuille law gives the flux $Q_{ij}$ in the channel at steady state:

\begin{equation} \label{eq:HP}
Q_{ij} = D_{ij}\frac{p_i - p_j}{L_{ij}}
\end{equation}

in which $p_i$ is the pressure at node $i$ and $D_{ij}$ is the conductivity of the channel $(i,j)$ given by:

\begin{equation} \label{eq:Dij}
D_{ij} = \frac{\pi r_{ij}^4}{8\eta}
\end{equation}

in which $\eta$ is the dynamic viscosity of the fluid. Note that $r_{ij} = r_{ji}$, and thus $D_{ij} = D_{ji}$.
If $Q_{ij} > 0$, then the fluid flows from $i$ to $j$; if $Q_{ij} < 0$, then the fluid flows in the opposite direction, from $j$ to $i$. As $D_{ij} = D_{ji}$ and $L_{ij} = L_{ji}$, then the direction of the flow (and the sign of $Q_{ij}$) is determined by the sign of $p_{ij} = p_i - p_j$.

H-P flow is the base of the mathematical formalism that will be used in this thesis (see section <ref>).

§.§ Adaptive Hagen-Poiseuille model derivation

The model used to describe the adaptive Hagen-Poiseuille flow on graphs will be the one described in <cit.>. This model consists on a class of equations that describe the flow of a viscous fluid through a network of straight channels (which can be described by a graph) with several static sources and sinks and adaptive channel conductivities.

This class of equations was derived as a way to model the dynamics and optimisation processes of Physarum polycephalum vein networks, and to correct previously derived models. Over the last twenty years, many other models were created aiming to describe Physarum's network's adaptive behaviour, namely vein formation and network optimisation over time. Some of these models also used the Hagen-Pouiseuille law to describe Physarum flow on graphs, and, to be able to achieve channel conductivity adaptation, they used local, specific evolution equations and stochastic update methods. The models converged to networks with steady paths, sometimes close to Steiner tree-like networks (see <cit.>) or shortest path solutions (see <cit.>). Despite their successes, these models were not consistent with some key physical notions: they did not preserve the volume of fluid of the network and, as such, the process of updating the conductivities on one site of the network was independent of the rest of the network. This is not a realistic mechanic as, for an incompressible fluid (which the H-P flow equation describes), for constant volume of fluid, any local change in conductivity must be followed by a global response. Thus, the class of equations used in this thesis preserves the volume of fluid and uses an update law that takes into consideration the entire network.

The class of equations used in this thesis and its derivation is now described. To describe a network of veins, one uses undirected graphs. Consider simple graphs $\mathcal{G}=(\mathcal{V},E)$ embedded in a $n$ dimensional Euclidean space. These graphs are made up of a set of $N$ nodes, $\mathcal{V}$, with spatial coordinates given by $ (x^1_i,...,x^n_N), \quad i = 1,...,N$, connected by a set of $M$ straight edges, $E$, described by $(i,j)$ (this describes the edge connecting nodes $i$ and $j$). The edges represent cylindrical, straight channels with elastic walls (that can expand or contract transversally to the flow of fluid). Thus, edge $(i,j)$ has a length $L_{ij}$ and radius $r_{ij}$. Each edge can also be described by its conductivity $D_{ij}$ (see eq. (<ref>)). The flow can be described using the Hagen-Poiseuille law (eq. (<ref>)).
The network has $K$ sources and $R$ sinks, which are located at fixed nodes. At these nodes, there is a inward flux of fluid to the network (in the case of a source) or an outward flux of fluid from the network (in the case of a sink). The inward/outward flux of fluid at node $j$ is represented by $S_j$. If $S_j > 0$, node $j$ is a source; if $S_j < 0$, node $j$ is a sink; if $S_j = 0$, node $j$ is a regular node.
As the fluid is incompressible,

\begin{equation} \label{eq:sources_sinks}
\sum_{j=1}^N S_j = \sum_{j:\text{sources}} S_j + \sum_{j:\text{sinks}} S_j = 0
\end{equation}

This shows the volume of fluid in the network should be constant. The volume of the cylindrical channel $(i,j)$ is given by $V_{ij} = \pi r_{ij}^2 L_{ij} = L_{ij}\sqrt{8\pi\eta}\sqrt{D_{ij}}$. As such, the (constant) volume of fluid in the network is given by:

\begin{equation} \label{eq:volume}
V = \sum_{(i,j)\in E} V_{ij} = \beta \sum_{(i,j)\in E} L_{ij}\sqrt{D_{ij}}
\end{equation}

with $\beta = \sqrt{8\pi\eta}$. Thus, the total volume of fluid $V$ is determined by the channel conductivities and the channel lengths. As the volume is constant, this volume is thus set by the initial conductivities of the system (see sec. <ref> for more information).
The steady state is determined by the Kirchhoff conservation laws, which state that flux is conserved at each node $j$:

\begin{equation} \label{eq:kirch}
\sum_{(i,j)\in E} Q_{ij} = S_i \quad,\quad i = 1, ... , N
\end{equation}

However, the fluxes and pressures determined by equations (<ref>) are a steady state solution obtained for fixed values of conductivities $D_{ij}$ and lengths $L_{ij}$. One wants to obtain an adaptive solution over time as, for living organisms, the flow of the veins is shown to alter the radii/conductivities of the channels <cit.>.

The adaptation law for the rate of change of the conductivities of the network $D_{ij}$ chosen was:

\begin{equation} \label{eq:Dij_ansatz}
\frac{d}{dt}\sqrt{D_{ij}} = f(\mathbf{Q}) - c\sqrt{D_{ij}}
\end{equation}

where $f(\mathbf{Q})$ is a function of the global flux of the network (with $f(\mathbf{0}) = 0$) and $c$ is some positive constant. Notice that $\sqrt{D_{ij}}$ is proportional to the transverse sectional area of channel $(i,j)$.

Using the condition of conservation of volume and eq. (<ref>), one gets:

$$\frac{d}{dt}V = \beta \sum_{(i,j)\in E}L_{ij}\frac{d}{dt}\sqrt{D_{ij}} = 0$$

and, using eq. (<ref>), one gets:

$$f(\mathbf{Q}) \sum_{(i,j)\in E}L_{ij} = c \sum_{(i,j)\in E} L_{ij} \sqrt{D_{ij}} = \frac{c}{\beta}V$$

which is constant.

One can define a new function $g$ by:

\begin{equation} \label{eq:f}
f(\mathbf{Q}) = V\frac{c}{\beta}\frac{g(Q_{ij})}{\sum_{(k,m)\in E} L_{km}g(Q_{km})}
\end{equation}

Substituting the function $f$ defined in (<ref>) in eq. (<ref>), one obtains:

\begin{equation} \label{eq:Dij_adapt}
\frac{d}{d\tau}\sqrt{D_{ij}} = \alpha\frac{g(Q_{ij})}{\sum_{(k,m)\in E}L_{km}g(Q_{km})} - \sqrt{D_{ij}} \quad,\quad (i,j)\in E
\end{equation}

where $\tau = ct$ and $\alpha = V/\beta = V/\sqrt{8\pi\eta}$.
For any choice of the function $g$, the adaptation equation (<ref>) ensures that the volume of fluid $V$ remains constant, and thus it contains a term that depends on the conductivities of all the edges of the network.

To be able to compute the temporal evolution of the system given by eq. (<ref>) and observe the dynamic changes in network shape and conductivity values, one needs to choose a $g$ function. This function was obtained by introducing a criterion based on minimizing dissipated energy, keeping the volume of fluid constant.
The dissipated power of the H-P steady flow $\mathcal{P}$ is given by:

\begin{equation} \label{eq:dissipated-power}
\mathcal{P} = \sum_{(i,j)\in E} \frac{Q_{ij}^2}{D_{ij}}L_{ij}
\end{equation}

Thus, using eqs. (<ref>) and (<ref>), one can construct the Lagrangian:

$$\mathcal{L} = \mathcal{P} - \lambda\left(V - \beta\sum_{(k,m)\in E}L_{km}\sqrt{D_{km}}\right)$$

where $\lambda$ is a Lagrange multiplier. Minimising $\mathcal{L}$ with respect to $D_{ij}$ and $\lambda$, and using eq. (<ref>), one gets:

\frac{\partial \mathcal{L}}{\partial \sqrt{D_{ij}}} & = \left(-\frac{Q^2_{ij}}{D^2_{ij}}L_{ij} + 2 \sum_{(k,m)\in E} \frac{Q_{km}}{D_{km}}\frac{\partial Q_{km}}{\partial D_{ij}}L_{km} + \lambda\beta\frac{L_{ij}}{2\sqrt{D_{ij}}}\right)2\sqrt{D_{ij}} = 0\\
\frac{\partial \mathcal{L}}{\partial \lambda} & = V - \beta\sum_{(k,m)\in E}L_{km}\sqrt{D_{km}} =
\end{dcases}$$

The second term on the right-hand side of the first equation is zero (compare to lemma 2.1 of <cit.>). As such, solving the equations in respect to $D_{ij}$ and $\lambda$, the conductivity values that minimise the dissipated power are:

\begin{equation} \label{eq:solution_lagrangian}
D^*_{ij} = \alpha^2\frac{Q^{4/3}_{ij}}{\sum_{(k,m)\in E}L_{km}Q^{2/3}_{km}} \qquad,\qquad D^{**}_{ij} = 0
\end{equation}

Comparing $D^*_{ij}$ of eq. (<ref>) to eq. (<ref>), one obtains the form of $g$ that minimises dissipated energy per unit time:

\begin{equation} \label{eq:g_23}
g(Q_{ij}) = Q_{ij}^{2/3}
\end{equation}

Thus, substituting eq. (<ref>) in eq. (<ref>) one obtains the adaptation equation used in this thesis:

\begin{equation} \label{eq:adaptation_eq}
\frac{d}{d\tau}\sqrt{D_{ij}} = \alpha\frac{|Q_{ij}|^{2/3}}{\sum_{(k,m)\in E}L_{km}|Q_{km}|^{2/3}} - \sqrt{D_{ij}} \quad,\quad (i,j)\in E
\end{equation}

where $\tau = ct$ and $\alpha = V/\beta$.

§.§ Murray's law

Murray's law is a theoretical law for an idealised network structure used in biophysical fluid dynamics that relates the radii of fluid-transporting, cylindrical veins at bifurcations in a network. It states that:

\begin{equation} \label{eq:murray_theo}
R^3 = \sum_i r_i^3
\end{equation}

in which, at a bifurcation in a network, $R$ is the radius of the parent branch and $r_i$ are the radii of the children branches.

Murray states that the total work involved in operating a section of a vessel should be minimum for an ideal network. This energy required to maintain the flow is given by the work required to overcome viscous drag forces (for a fluid described by the Hagen-Poiseuille law (<ref>), this work is given by $Q_{ij}\left(p_i - p_j\right)$) and the metabolic cost of maintaining the network structure, which is said to be proportional to its volume <cit.>.

Murray believed this idealised law would work on biological systems because living beings have undergone millions of years of Darwinian evolution, and have thus optimized their transport networks to an idealised state. Murray’s law has, in fact, been supported by measurements of small vessel branching relationships of arteries, in which the exponents obtained were in the range [2.7, 3.0] <cit.>, and by measurements of Physarum networks, where the measured exponents were in the range [2.53, 3.29] <cit.> (by "exponents", one means the exponents of equation (<ref>)).

The class of equations used in this thesis (described in section <ref>) uses similar energy minimisation considerations to the original Murray's law analysis. A significant difference between the two is that the formalism described takes the metabolic cost as a constraint (constant volume). Nevertheless,  a cubic law relationship between the flow rates and the radii of the channels was found ($Q_{ij} \propto D^{3/4}_{ij} \propto r^3_{ij}$, as seen in eq. (<ref>)). As the flux is conserved at each node (eq. (<ref>)), one obtains the generalised Murray's law <cit.>:

\begin{equation} \label{eq:generalised_murray}
\sum_{j:(i,j)\in E, Q_{ij} > 0} r_{ij}^3 = \sum_{j:(i,j)\in E, Q_{ij} < 0} r_{ij}^3 \quad,\quad i = 1, ... ,\text{N}
\end{equation}

Murray's law was verified for steady state flow using the class of equations described in sec. <ref> in <cit.>.

In this thesis, using the base formalism mentioned previously, I intend to verify if Murray's law is verified dynamically for adaptive flow with fixed sources and sinks.

§ STEINER TREE PROBLEM AND PHYSARUM

The Steiner minimum tree problem is a combinatorial optimisation problem. In its most general form, it asks for, given a set of objects, the optimal interconnect of all said objects according to a predefined objective function. When applied to undirected graphs, the Steiner tree problem seeks to find the tree that connects all nodes of interest (being able to contain nodes additional to the nodes of interest) and that also minimizes the total weight of the tree edges.

It's important to distinguish Steiner trees from minimum spanning trees (MST). For undirected graphs, the minimum spanning tree also seeks to find the tree that connects all nodes of interest and that also minimizes the total weight of the tree edges, but it cannot contain additional vertices besides the nodes of interest, unlike the Steiner minimum trees.

For graphs that symbolize paths between physical, geographical locations, like the graphs in this thesis, the "weight" that is to be minimized by the minimum Steiner tree is the length of said tree. Thus, for this thesis, the Steiner minimum tree (SMT) symbolizes the shortest tree that connects all nodes of interest.

Steiner minimum trees (in blue) for three regular polygons: a) a triangle, b) a square and c) a pentagon. Images courtesy of Martin Janecke, obtained from <cit.>.

As an example, figure <ref> shows the SMT (in blue) for three regular polygons. This figure also shows a key feature of SMTs: Steiner points, which are nodes that arise in the tree that are not part of the original group of nodes of interest. In figure <ref>, the Steiner points are the unmarked nodes that connect three (or more) edges of the blue SMT. For the triangle, there is a single Steiner point that corresponds to the center of said triangle; for the square, there are two Steiner points in the SMT and, for the pentagon, there are three Steiner points in the SMT.

Finding the Steiner minimum tree of a graph has a great deal of benefits and is a frequent problem in areas that involve building a cost-efficient and robust network of paths, like transportation, circuitry and telecommunications. After all, as a Steiner tree has a minimum length, it also has a minimum building and maintenance cost. Finding the SMT of a problem can also help increase the cost-efficiency of the network or can make it more robust (by providing minimum cost redundant paths).

Most versions of the Steiner problem in graphs are NP-hard <cit.> (this means it is a much harder problem than problems that can be solved by a nondeterministic Turing machine in polynomial time). Essentially, this means that it is unlikely that one is able to find an efficient algorithm to solve the Steiner tree problem for large graphs. As such, this is still an open problem. Biologically-inspired algorithms have thus tried to provide efficient ways to solve this very pertinent problem.

Real Physarum's network trees have already been shown to be able to converge to Steiner minimum trees, thus estimating solutions close to that of Steiner tree problems (figure <ref>) <cit.>. In figure <ref>, notice the resemblance between figures a3 and the SMT of the triangle; in figure <ref>, Physarum's network also shows a resemblance to the theoretical Steiner tree of the system. Physarum's networks aren't perfect: the paths aren't perfectly straight and the Steiner points aren't located at the most efficient locations. However, the Physarum networks shown in figure <ref> showed larger average values of FT/TL (fault tolerance / total length) than the the theoretical SMT values <cit.>. Fault tolerance is the probability of the network becoming disconnected when one of the links is removed, while the total length of the network is a measure of the cost of the network. Thus, as the Physarum network has a larger FT/TL value than the Steiner minimimum tree, it is a more robust and cheaper network than the SMT. An important note is that, while Physarum's network can converge to Steiner-like trees, they do not converge to a MST every time (see figures b-d of figure <ref>; figure b of fig. <ref> shows Physarum's network converging to the triangle's perimeter, for example).

Comparison between Physarum polycephalum networks and their theoretical Steiner minimum tree counterparts. a) For 3 food sources (FS) located at the edges of a plate, the Physarum specimen is made to cover the whole plate. The images show the network configuration after 0hr (a1), 6hrs (a2) and 36hrs (a3) of the start of the experiment. Images b-d show three different typical networks obtained. The graphic on the right shows the food sources locations as black circles and the SMT in green. b) For 6 FS, the image shows Physarum's network after some time had passed. The graphic above shows the food sources locations as black circles and the SMT in green. Images adapted from <cit.>.

As Physarum polycephalum's networks were shown to be able to converge to Steiner tree-like configurations, biologically-inspired algorithms were developed to mimic Physarum's behavior with the objective to reproduce these short networks. The Physarum network models described in section <ref> have been shown to converge to Steiner tree-like configurations, namely the hydrodynamics-based model Physarum solver <cit.>. Several hydrodynamics-based models have been used to solve real-world problems, like the minimal exposure problem in wireless sensor networks <cit.>, finding subnetworks for drug repositioning <cit.>, designing communication networks <cit.> and identifying elements of cancer-related signaling pathways <cit.>. The possible applications for these algorithms are vast and important.

CHAPTER: PROPERTIES OF ADAPTIVE HAGEN-POISEUILLE FLOWS

In order to determine the value and usefulness of a model, one must test it to see if it can accurately reproduce real physical results of what it is trying to portray. That is what we want to aim to do in this chapter of this thesis with the model described in section <ref>. First, we describe the algorithms used to implement said model (sec. <ref>). Then, we study the length distributions of several steady-state graphs obtained (sec. <ref>) and attempt to dynamically validate Murray's law as described in section <ref> (sec. <ref>). Next, we look at the properties of some simple networks (sec. <ref>). Finally, we compare the results obtained here with Physarum polycephalum's networks (sec. <ref>).

§ METHODS AND ALGORITHMS

As mentioned previously, the Physarum polycephalum network adaptation model used in this thesis is the one described in section <ref>. This model was implemented with the algorithms described below.

§.§ Initializing the simulation

First, the network graph $\mathcal{G}$ is initialized. $\mathcal{G}$ is a planar graph embedded in the two-dimensional Euclidean space. The graph is constructed as follows: first, a square lattice is created, and nodes are placed on the vertices of the lattice; the lattice is made to have side length $1$ and contains $N' \times N'$ nodes, $N'$ being the number of nodes on each side of the square. Then, the nodes' positions are randomly perturbed with Gaussian noise with a standard deviation $\sigma = 0.5$. Finally, the edges of the network result from a Delaunay triangulation between the nodes. Note that a planar Delaunay triangulation connects a set of nodes using triangular shapes in a way that no node is inside the circumscribed circle of any triangle of the triangulation.
The lengths of each edge $L_{ij}$ are calculated based on the positions of nodes $i$ and $j$ and are fixed for each different $\mathcal{G}$. In this thesis, only two different $\mathcal{G}$ of this kind will be taken into account, with $N' = \{25, 40\}$.

Some of the nodes of the network $j$ may be assigned as sources or sinks, having $S_j \neq 0$, and providing that condition of eq. (<ref>) is verified for the entire network.

The flow passing through the edges of the network $Q_{ij}$ is initialized by the starting conductivities $D_{ij}(t=0)$; these starting conductivities are initialized as positive random numbers with a homogeneous distribution.

Some parameters are kept constant for all runs, namely $\beta = 1$, $\Delta\tau = 0.1$, $V = 100$ and $\sum_{j\in\{\text{sources}\}}S_j = 1$ (unless stated otherwise).

§.§.§ Renormalization of volume

In order to be able to compare different runs of a simulation, some baseline conditions must be met for all the runs; namely, the volume of fluid should be the same for all the runs. While the model used in this thesis keeps the volume of the fluid in the network constant over time <cit.>, no measures were put in place to actually determine that the volume of fluid in each simulation was the same. Thus, we introduce a volume renormalization in this work to bridge this gap present in the previous description of the method.

As mentioned previously, for a given distribution of sources and sinks in a Delaunay graph, with prescribed input flows in the sources, the flow is initialized by the starting conductivities $D_{ij}(t=0)$. As eq. (<ref>) states that $V = \beta \sum_{(i,j)\in E} L_{ij}\sqrt{D_{ij}}$, the total volume of fluid in the graph is determined by $\beta$ (which is kept the same for all simulations), $L_{ij}$ (which is the same for simulations using the same baseline graph $\mathcal{G}$) and $D_{ij}$. Thus, the initial distribution of conductivities determines the total volume of fluid in the graph, and therefore the value of the parameter $\alpha$ in eq. (<ref>). The steady states obtained will depend on $\alpha$. Thus, for a different choice of initial conductivities, a different volume will be used and a different steady-state will be obtained.

As such, to analyze different steady states for the same volume of fluid, the initial conductivities must be renormalized. Consider some initial distribution of conductivities $\{D_{ij}(t=0) , (i,j)\in E\}$, which correspond to volume $V_0 = \beta \sum_{(i,j)\in E} L_{ij}\sqrt{D_{ij}}$. By multiplying each $D_{ij}(t=0)$ by a constant, that is, by rescaling the initial conductivities of the network, one can change the total volume of fluid in the network from $V_0$ to the desired volume $V$. The conductivities that correspond to volume $V$ are:

\begin{equation} \label{eq:renormalization_dij}
D_{ij}'(t=0) = \left(\frac{V}{V_0}\right)^2 D_{ij}(t=0)
\end{equation}

As such, to obtain the desired volume of fluid in the network $V$, one must, after initializing the conductivities with some random positive numbers, first calculate the current volume of fluid in the network $V_0$ and, afterward, multiply all $D_{ij}$ values by $(V/V_0)^2$.

With this renormalization, it is possible to compare different steady states for different initial conductivities and the same fluxes of sources and sinks.

§.§ Computing the temporal evolution of flows and conductivities

To describe the temporal evolution of the radii of the channels and the flow of fluid, two steps are repeated until a steady state is reached: first, the flows of the network are calculated by solving the linear equations described in eq. (<ref>), and then all the channel conductivities are adapted according to eq. (<ref>).

§.§.§ Calculation of the pressures and flows

Given a set of edge conductivities $D_{ij}$ and inward/outward fluxes of fluid at each node $S_j$, the linear equations described in eq. (<ref>) are solved to obtain the pressures at each node $p_i$ and the flows at each edge $Q_{ij}$. The method behind solving these equations will be the same as used in <cit.>. It is described below.

Substituting the H-P definition of flow of eq. (<ref>), one can rewrite each of the $i$ eqs. (<ref>) as:

\begin{equation} \label{eq:Cij}
\sum_{(i,j)\in E} \frac{D_{ij}}{L_{ij}}(p_i - p_j) = S_i \Leftrightarrow \sum_{(i,j)\in E} C_{ij}(p_i - p_j) = S_i
\end{equation}

in which $C_{ij} = D_{ij} / L_{ij}$ can be seen as weights of each of the edges of the network $(i,j)$ that measure the efficiency of each edge at carrying flow. One can re-write the left-hand side of eq. (<ref>) as:

\begin{equation}
\begin{aligned}
\sum_{(i,j)\in E} C_{ij}(p_i - p_j) & = \left(\sum_{(i,j)\in E} C_{ij}\right)p_i - \sum_{(i,j)\in E} C_{ij} p_j = \\
& = \left(\sum_{(i,k)\in E} C_{ik}\right)\sum_{(i,j)\in E}\delta_{ij}p_j - \sum_{(i,j)\in E} C_{ij} p_j = \\
& = \sum_{(i,j)\in E}\left[\left(\sum_{(i,j\in E)} C_{ik}\right)\delta_{ij} - C_{ij}\right]p_j = \\
& = \sum_{(i,j)\in E} w_{ij} p_j
\end{aligned}
\end{equation}

where $\delta_{ij}$ is the Kronecker delta. One can define a $N \times N$ symmetric matrix $\mathbf{W}$ with entries $w_{ij} = \left(\sum_{(i,j\in E)} C_{ik}\right)\delta_{ij} - C_{ij}$; this is the generalized Laplacian matrix of the network, as the network is a weighted graph with weights $C_{ij}$. As our network's graph is simple and undirected, this Laplacian matrix can be easily calculated by $\mathbf{W} = \mathbf{D} - \mathbf{A}$, where $\mathbf{D}$ is the degree matrix (a diagonal matrix whose $i$th diagonal entry is the sum of $C_{ij}$ for all the nodes $j$ connected to node $i$) and $\mathbf{A}$ is the adjacency matrix (a symmetric matrix with a null diagonal whose $ij$th entry is $C_{ij}$) of the network.
If $\mathbf{S}$ is the N-dimensional vector whose $i$th entry is $S_i$ (vector that describes the sources and sinks of the network) and $\mathbf{p}$ is the N-dimensional vector whose $i$th entry is $p_i$ (vector that describes the pressures of the nodes of the network), then one can rewrite the system of equations (<ref>) as:

\begin{equation}
\mathbf{W} \mathbf{p} = \mathbf{S}
\end{equation}

It's important to note that $\mathbf{W}$ is a singular matrix has all its rows sum up to zero, which means that it has an eigenvector with a null eigenvalue ($\mathbf{W}\mathbf{1} = \mathbf{0}$). Thus, $\text{rank}(\mathbf{W}) = N - 1$. This relates to the fact that, in the system of linear eqs. (<ref>), pressures are defined up to an additive constant (physically, one can only measure pressure differences). To be able to invert the matrix and then solve the system of equations, one adds a small arbitrary constant to one diagonal element of $\mathbf{W}$ (see section 2.4 of <cit.>). In the case of this thesis, each time equation (<ref>) is solved, one diagonal entry pertaining to a sink node is randomly perturbed. <cit.>

The pressure values are obtained by solving equation $\mathbf{p} = \mathbf{W}^{-1}\mathbf{S}$ using the Cholesky decomposition method, and the flow values are obtained by then using the H-P law (<ref>); these flow values are well-defined.

§.§.§ Dynamic adaptation of conductivities

After obtaining the $Q_{ij}$ flow values, all the channel conductivities are adapted according to eq. (<ref>).

In this thesis, the simple Euler method is used to numerically solve eq. (<ref>) using a time step $\Delta\tau = 0.1$. As such, the equation used to obtain the conductivities at time $t + \Delta\tau$ is:

\begin{equation} \label{eq:euler_method}
\sqrt{D_{ij}(t + \Delta\tau)} = \sqrt{D_{ij}(t)} + \Delta\tau\left(\frac{V}{\beta}\frac{\left|Q_{ij}(t)\right|^{2/3}}{\sum_{(k,m)\in E}L_{km}\left|Q_{km}(t)\right|^{2/3}} - \sqrt{D_{ij}(t)}\right) \quad,\quad (i,j)\in E
\end{equation}

where $Q_{ij}(t)$ denotes the flow value of edge $(i,j)$ at time $t$ (and similarly for $D_{ij}(t)$). Note that the time since the beginning of the simulation is given by $t = n\Delta\tau$, where $n$ is the number of iterations that occurred up to that point.

It's important to determine whether or not this method conserves the volume of fluid in the network over time. The volume of fluid at time $t + \Delta\tau$ is: (according to eq. (<ref>) and using eq. (<ref>))

\begin{equation}
\begin{aligned}
V(t+\Delta\tau) & = \beta\sum_{(i,j\in E)}L_{ij}\sqrt{D_{ij}(t+\Delta\tau)} = \\
& = \beta\sum_{(i,j\in E)}L_{ij}\sqrt{D_{ij}(t)} + \tau\left(V - \beta\sum_{(i,j)\in E}L_{ij}\sqrt{D_{ij}(t)}\right) = \\
& = V(t)
\end{aligned}
\end{equation}

Thus, equation (<ref>) preserves the volume of fluid of the network while also dynamically adapting the network's conductivity values.

§.§.§ Stopping criteria

For a static configuration of sources and sinks, the two steps described above (calculating the flows and then adapting the conductivities) are repeated over time until a steady state of the channel conductivities is reached. The numerical condition for reaching the steady state is:

\begin{equation} \label{eq:stop_condition}
\text{max}_{(i,j)\in E} \left|D_{ij}(t) - D_{ij}(t - \Delta\tau)\right| \le D_{\text{thresh}}
\end{equation}

where $t$ is the first time instance for which the condition is verified and $D_{\text{thresh}}$ symbolizes the desired precision value, being a positive small number. For this work, $D_{\text{thresh}} = 5\times10^{-6}$.

For non-static configurations of sources and sinks, the stopping criteria is different. By non-static configuration, one means that the sources and/or sinks locations change as time passes (at every time step of the simulation in the case of this work, unless stated otherwise). As the active sites change so frequently, the values of the conductivities $D_{ij}$ never get to a stable value, and change greatly every time step. However, for the cases considered in this thesis, the shape of the tree obtained does reach a constant form. As the topology of the network is what one wishes to study for these types of simulations, the stopping criteria was defined as follows: if the length of the network (defined by equation (<ref>)) remains unchanged for $N_{\text{iter}}$ iterations, then it is considered that the shape of the network was kept constant for that number of iterations and it is considered that the simulation reached its steady state. $N_{\text{iter}}$ is a parameter that is adjusted depending on the configuration of sources and sinks and the graph $\mathcal{G}$.

The simulations present in this work were carried out using Python, with the help of packages like  (for graph representation and analysis),  and  (for numerical computations) and  (for graphical representations and plots). The code used was originally developed by Almeida <cit.> and was adapted for this work.

§ LENGTH DISTRIBUTION OF THE STEADY STATE GRAPHS

As mentioned previously (in sec. <ref>), different choices of $\{D_{ij}(t=0) , (i,j)\in E\}$ may lead to trees with different shapes connecting sources and sinks. These trees are steady states of the H-P flow.

In some cases, we may have several apparently disconnected trees. However, the graph tree has only one connected component at all times, with conductivities partitioned in two sets of high and low values. This allows for one to identify the channels that are effectively conducing flow in the tree. <cit.>

To be able to evaluate, characterize and quantify the different possible steady states, different parameters are used. The first parameter studied will be the length of the tree, measured at steady state. The total length of the graph is given by:

\begin{equation} \label{eq:total_length}
L = \sum_{(i,j)\in E'} L_{ij}
\end{equation}

where $E'$ is the set of effectively conducting edges, which are edges with conductivity values above a certain threshold, that is, $E' = \{(i,j) \in E: D_{ij} > D_{\text{thr}}\}$. In the case of this work, $D_{\text{thr}} = 5 \times 10^{-4}$.

$25\times25$ lattice. The mean of the length distribution is $\overline{L} = 1.311$ and the standard deviation is $\sigma = 0.030$.

$40\times40$ lattice. The mean of the length distribution is $\overline{L} = 1.247$ and the standard deviation is $\sigma = 0.017$.

Length distribution of the steady state trees (right), shortest steady state tree obtained (left) and longest steady state tree obtained (middle), for one source and one sink, for different $\mathcal{G}$. The length of the steady state tree is written above its image. The linear distance between source and sink is $L_{\text{lin}} = 4\sqrt{2}/5 \approx 1.131$. In principle, in the limit $N' \to \infty$, $L\to L_{\text{lin}}$.

A study of the length distribution of the steady state graphs for different sets of sources and sinks can be seen in figures <ref> and <ref> for one source and one sink and five sources and four sinks, respectively. The trees were generated on the square lattice graph described previously of side length $1$ with $N = N'\times N'$ nodes, with $V = 100$ and $\sum_{j\in\text{sources}}S_j = 1$. Each figure contains the results for a graph $\mathcal{G}$ with $N' = 25$ and $N' = 40$. For each set of sources and sinks and for each $\mathcal{G}$, 300 trees were generated, each with different random initial conductivity conditions.

For all the steady state figures presented in this paper, the initial Delaunay triangulation $\mathcal{G}$ is shown in grey, and the thickness of the black line is proportional to the radius of the edges $\sim D_{ij}^{1/4}$. Additionally, source nodes are represented as circles, and sink nodes are represented as triangles.

The length distributions of both figures <ref> and <ref> resemble normal distributions. These images confirm the statement that different initial conductivity conditions lead to several different steady states.

As the size of the lattice $N'$ increases, the continuous limit is approached. For the case of one single source and one single sink (figure <ref>), the shortest graph connecting them tends to reach a straight line between the source and the sink (compare figs. <ref> to <ref>). This is corroborated by the fact that the length of the shortest tree for the $N' = 40$ lattice is smaller than the length of the shortest tree for $N' = 25$, while also being closer to the linear distance between source and sink ($L_{25} > L_{40} > L_{\text{lin}}$). This is also true for the mean value of the length distribution of both lattices ($\overline{L_{25}} > \overline{L_{40}} > L_{\text{lin}}$).

The tendency of the simulation to take the shortest path possible between source and sink can also be shown by the fact that the standard deviation of the length distribution of the $N' = 40$ case is smaller than that of the $N' = 25$ case. This is because, when the number of nodes increases, so does the number of edges and therefore the number of possible paths the fluid can take; a smaller standard deviation shows the fluid didn't create very long trees and stayed as close to the shortest path as possible.

$25\times25$ lattice. The mean of the length distribution is $\overline{L} = 3.40$ and the standard deviation is $\sigma = 0.27$.

$40\times40$ lattice. The mean of the length distribution is $\overline{L} = 3.44$ and the standard deviation is $\sigma = 0.33$.

Length distribution of the steady state trees (right), shortest steady state tree obtained (left) and longest steady state tree obtained (middle), for five sources and four sinks, for different lattices. The length of the steady state tree is written above its image.

For the case of five sources and four sinks (figure <ref>), we once again observe the tendency of larger node number leading to shorter tree lengths, as once again the length of the shortest tree for the $N' = 40$ lattice is smaller than the length of the shortest tree for $N' = 25$. However, this time both the mean and the standard deviation of the length distribution for the $N' = 40$ lattice are larger than that of the $N' = 25$ case. This might be because of the fact that with more sources and sinks come a lot more possibilities of topologies of the tree, that can be further expressed with a larger node number.

The case of two sources and two sinks was also studied and, to illustrate it, two different trees were generated and are presented in figure <ref>. Even though the source and sink nodes were the same, each tree was generated with different values of incoming/outgoing flow of fluid at each source/sink node $j$, $S_j$. The trees were generated using the same initial conditions, and were generated on the square lattice graph described previously of side length $1$ with $N = 25\times 25$ nodes, with $V = 100$ and $\sum_{j\in\text{sources}}S_j = 1$.

$S_{\text{left source}} = S_{\text{right source}} = 0.5 = - S_{\text{left sink}} = - S_{\text{right sink}}$.

$S_{\text{left source}} = - S_{\text{right sink}} = 0.3$, and $S_{\text{right source}} = - S_{\text{left sink}} = 0.7$.

Two different possible steady state configurations for two sources and two sinks in a $25\times25$ lattice.

Figure <ref> shows that not all steady state trees are necessarily visibly connected. Figure <ref> seems like it contains two disconnected graphs. However, as stated earlier, the graph tree has only one connected component at all times, with conductivities partitioned in two sets of high and low values, which allows the identification of the channels that are effectively conducing flow in the tree.

§ VALIDATION OF MURRAY'S LAW

While the generalised Murray's law (eq. (<ref>)) was verified for steady state mathematically, it has yet to be verified in simulation. As it is a direct consequence of the set of equations used to determine the steady state, the hypothesis that it is verified at steady state is likely, but it is unknown whether or not it is verified dynamically.

To determine if the generalised Murray's law (<ref>) is verified, the configuration of sources and sinks of figure <ref> will be used. First, the initial conductivities will be chosen to be $D_{ij}(t=0) = 1, \forall (i,j)\in E$, in order to ensure the steady state obtained will always be the same for all the simulations. Then, the steady state will be observed and two separate bifurcation points in the tree (at steady state) will be identified. By bifurcation points, one means nodes $j$ which are connected to more than two conducting edges. For these two chosen bifurcation nodes, the quantity $\left|\sum_{j:(i,j)\in E, Q_{ij} < 0} r_{ij}^3 - \sum_{j:(i,j)\in E, Q_{ij} > 0} r_{ij}^3\right| := \left|R_{\text{out}} - R_{\text{in}}\right|$ will be measured as time progresses. For Murray's law to be verified, the quantity $\left|R_{\text{out}} - R_{\text{in}}\right|$ must be zero.

The graph of $\left|R_{\text{out}} - R_{\text{in}}\right|$ (normalized to its initial value) as a function of time is shown for both bifurcation nodes chosen in figure <ref>, as well as the steady state and the bifurcation nodes chosen (the bifurcation nodes chosen are surrounded by a black square box).

Steady state analyzed, with the bifurcation nodes chosen surrounded by black square boxes.

Variation of $\left|R_{\text{out}} - R_{\text{in}}\right|$ over time for leftmost bifurcation node.

Variation of $\left|R_{\text{out}} - R_{\text{in}}\right|$ over time for rightmost bifurcation node.

Dynamic verification of Murray's law for two different bifurcation nodes.

As the graphs of figures <ref> and <ref> show, $\left|R_{\text{out}} - R_{\text{in}}\right| > 0$ for some initial portion of time for both nodes; after some time has passed, at steady state, $\left|R_{\text{out}} - R_{\text{in}}\right| = 0$. Thus, the generalized Murray's law (<ref>) is not verified dynamically but is verified at steady state.

§ PROPERTIES OF SIMPLE NETWORKS

§ COMPARISON WITH PHYSARUM POLYCEPHALUM

CHAPTER: PATHS OF EXTREMAL LENGTHS CONNECTING POINTS IN AN EUCLIDEAN SPACE

The goal of this chapter is to analyze extremal length graphs connecting $n$ points in an Euclidean space, using the properties of adaptive Hagen-Poiseuille flows. We will explore the ability of the model studied in this work to generate approximate Steiner-type solutions in graph theory. We will develop several strategies for the design of transportation systems and we propose an algorithm to determine their efficiency.

§ FINDING PATHS OF EXTREMAL LENGTHS CONNECTING THE NODES OF A GRAPH

The survival of an organism is greatly dependent on its biological network. For optimal survivability, a biological network must be efficient in transporting and distributing resources and nutrients throughout the body (by providing a short path between relevant sites). In addition to ensuring transport efficiency, a biological network should also be as short in length as possible, so as to minimise the costs of building said network. Furthermore, biological networks should be robust and redundant as well, to prevent network failure in the case of damage to some channel. In the case of Physarum, its network's adaptability is a key factor to its survivability.

These factors are not only relevant in biological networks; in fact, they are key concepts behind the design of many man-made networks, transport networks (railroads, highways) being an example.

To first illustrate these concepts, we will present geometric examples. We consider a set of $n$ points in an Euclidean space of dimension $\mathcal{N}$. These points may be connected through a connected graph, and there is an infinity of solutions for the connected graph.

Let us first consider the case of $n = 3$ points in an equilateral triangular configuration on a 2D plane. Two possible trees for the graph connecting these three points are shown in figure <ref>.

Minimal Steiner tree of a triangle.

Perimeter of a triangle.

Two possible connections of 3 points in a 2D graph.

The points $p_i, i=\{1, 2, 3\}$ make an equilateral triangle of side 1, while the point $p_S$ is the center of the same triangle. We assume the following coordinates for the 3 vertices and central point:

$$p_1 = (0,0) \qquad p_2 = (1,0) \qquad p_3 = \left(\frac{1}{2},\frac{\sqrt{3}}{2}\right) \qquad p_S = \left(\frac{1}{2},\frac{\sqrt{3}}{4}\right)$$

The length of the tree of figure <ref> is $L_a = 3\sqrt{7}/4\approx1.98$, while the graph length of the graph of figure <ref> is $L_b = 3$. The graph of figure <ref> is the graph of minimal lengths that connects the three points (and that is allowed to contain additional points, $p_S$), and thus is the Steiner minimum tree of the triangle vertices.

While this is an abstract system, the tree connecting the three points can also symbolize a multitude of physical, social or biological systems; for example, it can be a model of a system of roads connecting three cities. Suppose each road is used by frequent travellers that go from city $i$ to city $j$, and that the number of travelers that go from $i$ to $j$, for every $(i,j)$ pair, is asymptotically the same. For the road system described by the graph of figure <ref>, the mean distance between two cities is $M_a = 1.98$, and for the road system described by the graph of figure <ref>, this distance is $M_b = 1$.

Therefore, from the point of view of communications, the graph configuration of figure <ref> is more advantageous in the sense that is has lower transportation costs for the travelers. However, the graph configuration of figure <ref> is more advantageous in the sense that it has lower building and maintenance costs.

For a given set of configurations of points in an Euclidean space, we want to decide which graph fulfills the requirement of efficiency and low costs.

§.§ Stochastic algorithms used

As mentioned in <ref>, different initial conductivity values will lead to different steady-state tree shapes. To obtain different trees with different shapes, efficiencies, and lengths, a stochastic search will be done. The simulation will be run $N_{\text{runs}} = 300$ times, each time with different random initial conditions, for each configuration with $n$ sites.

The adaptive Hagen-Poiseuille flows equations require that there is at least one source node and one sink node for a tree to form. In this paper, for the case of a configuration with $n$ sites, each site is equal to every other site, making each site equally likely to be a source or a sink. As such, different stochastic algorithms will be used to choose which nodes are sources and which nodes are sink at each iteration step.

Each algorithm is given a name for quick identification. The algorithms are the following: at each iteration (every $\Delta t$), let $I_0 = \sum_{i\in\text{sources}} S_i$ be the total inward/outward flux of fluid to/from the network:

* : 1 source (with $S_i = I_0$) and 1 sink (with $S_j = -I_0$) are chosen randomly among the $n$ sites;

* : $n/2$ sources (each with $S_i = I_0/(\frac{n}{2})$) and $n/2$ sinks (each with $S_j = -I_0/(\frac{n}{2})$) are chosen randomly from the $n$ sites; if $n$ is an odd number, one of the sites is inactive;

* : 1 source (with $S_i = I_0$) is chosen randomly from the $n$ sites and the remaining $n - 1$ sites are sinks (each with $S_j = - I_0 / (n-1)$);

* : $1 \leq n_{so} \leq n - 1$ sources are chosen randomly from the $n$ sites (each with $S_i = I_0/n_{so}$), and the remaining $n-n_{so}$ sites are sinks (each with $S_j = -I_0/(n-n_{so})$).

§.§ Characterization parameters

Three metrics will be used to characterize each steady state: its length $L$, its average length between pairs of sites $\nu$, and its cost-efficiency CE. These three parameters will be measured for all $N$ steady-state trees.

\begin{equation} \label{eq:L}
L = \sum_{(i,j)\in E'} L_{ij}
\end{equation}

be the length $L$ of a graph, where $E'$ is the set of effectively conducting edges, that is, edges with conductivity values above a certain threshold $(E' = \{(i,j)\in E: D_{ij} > D_{\text{thresh}})\}$, with $D_{\text{thresh}} = 5\times10^{-4}$). Note that $L_{ij}$ is the length of edge $(i,j)$.

\begin{equation} \label{eq:nu}
\nu = \frac{1}{\left(^n_2\right)} \sum_{\substack{i,j\in\{\text{sites}\},\\i\neq j,\\ (i,j)\neq(j,i)}} L_{s_i s_j} = \frac{2}{n(n-1)}\sum_{\substack{i,j\in\{\text{sites}\},\\i\neq j,\\ (i,j)\neq(j,i)}} L_{s_i s_j}
\end{equation}

be the average length between pairs of sites $\nu$, where $L_{s_i s_j}$ is the length of the path between sites $s_i$ and $s_j$. The quantity $L_{s_i s_j}$ was computed by finding the shortest path between two sites and adding up the lengths of all conducting edges on said path; the shortest path between the two sites was found using an algorithm very similar to Breadth first search <cit.>.

\begin{equation} \label{eq:CE}
\text{CE} = \frac{1}{L\nu}
\end{equation}

be the cost-efficiency CE. By thinking once again of these steady-state trees as roads connecting cities, an efficient tree in terms of communication and travel would have a low $\nu$ value, while a tree with a low $L$ value would be cheap to construct or maintain. The parameter CE tries to take into account both of these facts.

§ OPTIMISING THE STEINER TREE SOLUTION

In section <ref>, the difference between the Steiner minimum tree and the minimum spanning tree has been explained: the STM is able to contain additional vertices other than the nodes of interest, called Steiner points, unlike the minimum spanning tree.

One of the characteristics of the steady-state solutions of the H-P adaptive flows is the existence of additional vertices other than sources and sinks (as discussed in section <ref> and seen in figures <ref>, <ref> and <ref>). In some instances, the steady state of the H-P adaptive flow has the topological structure of a Steiner tree, despite not needing to necessarily be a minimum-length tree. In this case, when the topology of the graph is fixed, it is possible to adjust the coordinates of the Steiner points in such a way that the steady-state graph becomes a Steiner minimum tree.

Trees that resemble Steiner minimum trees for three polygons: a) $n=3$, b) $n=4$ and c) $n=5$.

We exemplify this with simple examples. Consider the graphs in figure <ref>. All the trees have the topology of a Steiner tree, that is, they look like they could be a Steiner minimum tree, but their Steiner points may be slightly out of place.

Let's look at the case of figure <ref>, the triangle-like configuration. Consider the Steiner point is located at coordinates $(x_{S}, y_{S})$. The sum of the squares of the lengths of all the edges of the tree is

\begin{equation}
L^2 = \sum_{i=1}^3 (x_i - x_S)^2 + (y_i - y_S)^2,
\end{equation}

Minimizing $L^2$ relative to $x_S$ and $y_S$,

\begin{equation} \label{eq:minimizing_L_steiner}
\begin{aligned}
\frac{\partial L^2}{\partial x_S} & = - 2 \sum_{i=1}^3 (x_i - x_S) = 0 \\
\frac{\partial L^2}{\partial y_S} & = - 2 \sum_{i=1}^3 (y_i - y_S) = 0, \\
\end{aligned}
\end{equation}

we obtain the solution

\begin{equation} \label{eq:steiner_points_triangle_minimum}
x_S = \frac{1}{3} \sum_{i=1}^3 x_i \quad,\quad y_S = \frac{1}{3} \sum_{i=1}^3 y_i
\end{equation}

It's important to note that, since

\begin{equation}
\frac{\partial^2 L^2}{\partial x_S^2} = \frac{\partial^2 L^2}{\partial y_S^2} = 2 > 0
\end{equation}

then the solution described in eq. (<ref>) is in fact a minimum. As such, eq. (<ref>) gives the positions of the Steiner point for the tree of figure <ref>. In fact, eq. (<ref>) provides the positions of the Steiner point for all trees that possess that general shape, one Steiner point and that result from configurations of three sites.

This procedure is trivially generalised for any tree shape, number of sites and any number of Steiner points. All one must do is construct $L^2$ appropriately, taking into account all paths of the tree (paths that connect sites to Steiner points and paths that connect Steiner points themselves), and then obtain its minimum by calculating $\partial L / \partial x_{Sj}$ for all $j$ Steiner points and equating it to zero.

For example, for the tree of picture <ref> (the square-like configuration with 2 Steiner points), if the Steiner point $i$ has coordinates $(x_{Si},y_{Si})$ and the node of interest $j$ has coordinates $(x_j, y_j)$, $L^2$ is:

\begin{equation}
L^2 = \sum_{i=1}^4 \sum_{j=1}^2 \left((x_i - x_{Sj})^2 + (y_i - y_{Sj})^2\right) + (x_{S1} - x_{S2})^2 + (y_{S1} - y_{S2})^2
\end{equation}

After minimization, we obtain the new positions of the Steiner points:

\begin{equation}
\begin{aligned}
x_{S1} & = \frac{1}{8}(x_1 + x_2 + 3x_3 + 3x_4) \quad,\quad y_{S1} & = \frac{1}{8}(y_1 + y_2 + 3y_3 + 3y_4) \\
x_{S2} & = \frac{1}{8}(3x_1 + 3x_2 + x_3 + x_4) \quad,\quad y_{S2} & = \frac{1}{8}(3y_1 + 3y_2 + y_3 + y_4)
\end{aligned}
\end{equation}

If $x_1 + x_2 = x_3 + x_4$ and $y_1 + y_2 = y_3 + y_4$, the two Steiner points coalesce into a single Steiner point.

For the tree of picture <ref> (the pentagonal-like configuration with 3 Steiner points), the result obtained is:

\begin{equation}
\begin{aligned}
x_{S1} = \frac{1}{21}(8x_1 + 3x_2 + x_3 + x_4 + 8x_5) & \quad,\quad \text{and similarly for } y_{S1} \\
x_{S2} = \frac{1}{7}(x_1 + 3x_2 + x_3 + x_4 + x_5) & \quad,\quad \text{and similarly for } y_{S2}\\
x_{S3} = \frac{1}{21}(x_1 + 3x_2 + 8x_3 + 8x_4 + x_5) & \quad,\quad \text{and similarly for } y_{S3}\\
\end{aligned}
\end{equation}

Many times, the problem in calculating Steiner trees resides in the fact that it is unclear how many Steiner points there may be and how they could be connected to each other and to the sites. Using the algorithm of the model studied in this work, one can determine how many Steiner points the minimum tree should have (as the algorithms produce trees that tend to the shortest paths); then, using these calculations, one can determine the Steiner points' exact locations, based on the sites used and, in this case, this is a solvable linear problem.

§ RESULTS

The algorithms and parameters described will first be tested and studied on networks connecting vertices of regular polygons on a planar graph embedded in the two-dimensional Euclidean space. Then, the algorithms and parameters will be used to optimize paths for communication and transport systems; as an example, we present mainland Portugal.

§.§ Configurations of extremal lengths

To study the algorithms and parameters described, the methods described in section <ref> will be used to initialize the simulation and to compute the temporal evolution of flow and conductivities. Among the $N' \times N'$ nodes that make up the mesh, $n$ sites will be chosen as possible sources/sinks for the simulations.

The parameters used were: $\beta = 1$, $\sum_{i\in\text{active sources}}S_i = 1$, $V = 100$, $N' = 25$, $\Delta\tau = 0.1$ and $N_{\text{runs}} = 300$.

The cases of $n = \{3, 4, 5\}$ will be studied, with the sites disposed in a regular polygon configuration. The results are presented in the following sections. For additional results, see appendix <ref>.

§.§.§ Triangular configuration

For $n = 3$, the sites approximately form an equilateral triangle of side $0.8$. The theoretical value for the perimeter of the triangle is $P_\triangle = 2.275$ (calculated using the coordinates of the sites) and for the minimum Steiner tree length it is $L_{\text{Stei}_\triangle} = 1.313$ (calculated using the Steiner point location obtained with the equations of section <ref>).

The results for algorithms ,  and  are shown in figures <ref>, <ref> and <ref> respectively. Note that algorithms  and  are equivalent for the case of $n = 3$. The results shown are the steady states obtained with smallest $L$ and $\nu$ values and with the biggest CE value, and the graphs of $\nu$ as a function of $L$ for all runs and for all algorithms can be seen in figure <ref>. Additional results can be viewed in appendix <ref>.

The tree obtained as the minimum length, minimum $\nu$ and maximum CE tree for algorithms  and  is very similar to the minimum Steiner tree of the triangle. For these two algorithms, the graph of $\nu$ vs. $L$ shows a linear relation between the two variables: a tree with a large length would also present a large $\nu$, thus a lower efficiency.

For the algorithm , the minimum Steiner tree was not obtained; the minimum $\nu$ tree was the triangle with its full perimeter and the minimum $L$ and maximum CE tree was two sides of the triangle. For this algorithm, the graph of $\nu$ vs. $L$ shows two clusters: one representing trees with large $\nu$ and small $L$ (in which the two variables seem to show a linear relation) and one representing trees with small $\nu$ and large $L$ (where the lower bound of this region appears to present a linear relation).

The algorithm that produced the shortest tree and the tree with the largest CE value was , which is also the algorithm that most frequently produced short and high CE value trees (see fig. <ref> and table <ref>).

The algorithm that produced the smallest $\nu$ tree was , which was also the algorithm that most frequently produced small $\nu$ value trees (see fig. <ref> and table <ref>).

§.§.§ Square configuration

For $n = 4$, the sites approximately form a square of side $0.7$. The theoretical value for the perimeter of the square is $P_\square = 2.715$ (calculated using the coordinates of the sites) and, for the minimum Steiner tree length, it is $L_{\text{Stei}_\square} = 1.851$ (calculated using the Steiner points locations obtained with the equations of section <ref>).

The results for algorithms , ,  and  are shown in figures <ref>, <ref>, <ref> and <ref> respectively. The results shown are the steady states obtained with the smallest $L$ and $\nu$ and largest CE values, and the graph of $\nu$ as a function of $L$ for all runs for all algorithms can be seen in figure <ref>. Additional results can be viewed in appendix <ref>.

Algorithm  once again did not obtain the minimum Steiner tree for the square. The tree shape obtained was the full perimeter of the square. Once again, the variables $\nu$ and $L$ show a somewhat linear relation, but graph <ref> now displays points located in a region rather than showing a straight line of points.

The remaining algorithms, ,  and , all obtained a tree very similar to the minimum Steiner tree of the square for the trees with minimum length and maximum CE. The trees with minimum $\nu$ were the (slightly distorted) full perimeter of the square. Once again, the graph of $\nu$ vs. $L$ shows two clusters: one representing trees with large $\nu$ and small $L$ and one representing trees with small $\nu$ and large $L$. Very few runs fall between these two regions (less than 5% for all these algorithms).

It's relevant to note that the shortest tree is usually the tree with largest cost-efficiency value.

The algorithm that produced the shortest tree and the tree with the largest CE value was .
The algorithm that most frequently produced short and high CE value trees was  (see fig. <ref> and table <ref>).

The algorithm that produced the smallest $\nu$ tree was . The algorithm that most frequently produced small $\nu$ value trees was  (see fig. <ref> and table <ref>).

These solutions could potentially be improved by using the algorithm described in section <ref> to calculate the best position for the Steiner points.

§.§.§ Pentagonal configuration

For $n = 5$, the sites approximately form a pentagon of side $0.55$. The theoretical value for the perimeter of the pentagon is $P_{\text{penta}} = 2.802$ (calculated using the coordinates of the sites) and for the minimum Steiner tree length it is $L_{\text{Stei}_{\text{penta}}} = 2.194$ (calculated using the Steiner points locations obtained with the equations of section <ref>).

The results for algorithms , ,  and  are shown in figures <ref>, <ref>, <ref> and <ref> respectively. The results shown are the steady states obtained with smallest $L$ and $\nu$ values and largest CE value, and the graph of $\nu$ as a function of $L$ for all runs for all algorithms can be see in figure <ref>. Additional results can be viewed in appendix <ref>.

Algorithm  once again did not obtain the minimum Steiner tree for the pentagon. The result obtained was the full perimeter of the pentagon. Once more, the variables $\nu$ and $L$ show a linear relation.

Algorithms  and  obtained a tree very similar to the minimum Steiner tree of the pentagon for the tree with minimum length for the former and for the tree with maximum CE for the latter, while algorithm  obtained a tree which looks like a Steiner tree (right side) conjoined with the perimeter of the pentagon (left side) for the shortest length tree.
The tree with minimum $\nu$ for algorithm  was the (distorted) full perimeter of the square, while for algorithms  and  it was a mix of the Steiner tree of the pentagon with the perimeter of the pentagon. Once again, the graph of $\nu$ vs. $L$ shows two clusters for all these algorithms: one representing trees with large $\nu$ and small $L$ and one representing trees with small $\nu$ and large $L$.

The algorithm that produced the shortest tree was . The algorithm that most frequently produced short trees was  as well (see fig. <ref> and table <ref>).

The algorithm that produced the smallest $\nu$ tree was . The algorithm that most frequently produced small $\nu$ value trees was  (see fig. <ref> and table <ref>).

The algorithm that produced the tree with the largest CE value was . The algorithm that most frequently produced high CE value trees was  (see fig. <ref> and table <ref>).

§.§.§ Remarks about polygonal studies

Using the adaptive H-P model studied in this thesis, the algorithms studied in section <ref> were all able to obtain Steiner minimal trees for all polygonal configurations, with the exception of the  algorithm that was unable to produce a Steiner-like tree.

The tree shape obtained that corresponded to the smallest $\nu$ value, thus the highest efficiency, was mostly the perimeter of the polygon (often including a Steiner point near some of the sites). The tree shape obtained that corresponded to the smallest $L$ value, thus the lowest cost, was mostly the Steiner minimal tree-like result. The tree shape obtained that corresponded to the largest CE value was usually the same as the smallest $L$ value one, that is, the Steiner minimal tree-like result.

The algorithm that produced the shortest trees (smallest $L$ value) was consistently  for all configurations.
The algorithms that produced the most efficient trees (smallest $\nu$ value) were  for $n = \{3,5\}$ and  for $n = \{4\}$.
The algorithms that produced the trees with the highest cost-efficiency values were  for $n = \{3,4\}$ and  for $n = \{5\}$.

The algorithms that produced Steiner tree-like shapes more frequently were the  algorithm (for the square configuration) and the  algorithm (for the pentagonal configuration).
For the triangular configuration, algorithms  and  were equally likely to produce Steiner-like trees.
This is observed by analyzing the length distributions present in appendix <ref>.

Physarum, as seen in figure <ref>, is able to produce Steiner minimal tree-like networks, even though it does not produce one every single time. This is exactly what was observed for the simulations of section <ref>: while the minimal length steady state trees were obtained a few times for each configuration, they were not obtained every time. This is because the steady state tree shape is heavily dependent on the initial conditions of the conductivities $D_{ij}(t=0)$, as first seen in several simulations of chapter <ref>. Similarly, Physarum's networks are heavily dependent on external stimuli (like food sources and light sources, for example). Figure <ref> shows a Physarum specimen growing on an agar-covered surface; many external stimuli could lead Physarum to construct different shapes of its network, like non-uniformities in the agar, for example.

§.§ Optimizing paths for communication systems

As an example, the algorithms described will now be used to approximate the mainland Portuguese railroad transport network. This is a continuation of the work started in <cit.>, and thus we will use many of the same resources (like the original mesh, railway tree and MST) as the previous work.

The mesh used was first described in <cit.>. It is shaped like mainland Portugal and is made up of 1005 nodes and 2817 edges resultant of triangulation. We chose 25 sites of interest, which are the 18 Portuguese district capitals and 7 other major cities. These sites and the real railway tree that connects them can be seen in figure <ref>, as well as the MST that connects all sites. The length of each edge $L_{ij}$ used was the geodesic distance between nodes $i$ and $j$, which was calculated approximately using the Haversine formula.

Railway with stations.

a) Approximation of part of the railway of mainland Portugal that shows all 18 district capitals (in green) (note that Viseu's district capital was replaced by Mangualde) and 5 additional cities (in red). Some rail lines are no longer active but are included here for consistency. b) Tree that represents Portugal's railway system. c) Minimum spanning tree that connects all railway sites. Images adapted from <cit.>.

The parameters used were: $\beta = 1$, $\sum_{i\in\text{active sources}}S_i = 1$, $V = 100$, $\Delta \tau = 0.1$ and $N_{\text{runs}} = 300$.

The results for algorithms , ,  and  are shown in figures <ref>, <ref>, <ref> and <ref> respectively. The results shown are the steady states obtained with smallest $L$, $\nu$ and CE values, and the graph of steady state parameters $\nu$ as a function of $L$ for all runs (figure <ref>). Additional results can be seen in appendix <ref>.

$\nu$ as a function of $L$ for Portugal's configuration, for all algorithms and all runs. Legend: $\blacktriangle$: ; $\bullet$: ; $\Diamond$: ; $\blacksquare$: .

Before any analysis of the results takes place, it is important to refer that algorithm  was not like the others when it comes to running the simulations. It took about 10 times more to run 300 iterations than the other algorithms and only $67\%$ of iterations actually converged (using the stopping criteria of the length of the tree being unchanged for 500 iterations) with a maximum number of iterations of 5000. For this study, only the converged runs were taken into account. This makes  a much less efficient algorithm when compared to the others on a computational level, especially when the number of sites used increases.

First, let us look at the tendencies of each algorithm. Using figure <ref> and the images and table of appendix <ref>, we can see that algorithms  and  produce very similar results: the trees show a relatively concentrated cluster of runs with low $L$ and (mostly) low $\nu$. We can see that these algorithms show very concentrated $L$ values and much more dispersed $\nu$ values.
On the other hand, algorithm  shows a much more dispersed cluster of trees at approximately the same location (with low $L$) but reaching higher $\nu$ values. It shows the most dispersed $\nu$ values of all algorithms and $L$ values slightly more dispersed than the previously mentioned algorithms.
Algorithm  is the most disparate of them all, with a dispersed cluster of trees with high $L$ values and low $\nu$ values. While this cluster is dispersed in $L$ values (unlike the others), it has the lowest standard deviation with regard to $\nu$ values of all the algorithms.

Let us now comment on the shape of the best trees obtained.  shows the trees with the most numerous loops. The smallest $L$ tree, smallest $\nu$ tree, and largest CE tree produced by this algorithm present five, eight, and five loops, respectively. Compared with the other algorithms,  presents the fewest loops, with zero, three, and one loop,  presents two, five, and three loops, and  presents two, four and three loops, respectively. The real railway tree shows four loops, while the MST shows zero loops (see figure <ref>). We can see that the most efficient trees present the highest number of loops and the shortest trees the lowest number of loops, with the high CE trees showing a middle ground between the two.

The shortest tree obtained was produced using algorithm . This is also the algorithm that most frequently produces short trees (see fig. <ref> and table <ref>).

The tree with lowest $\nu$ value obtained was produced with algorithm . The algorithm that most consistently produces low $\nu$ trees is  (see fig. <ref> and table <ref>).

The tree with the highest CE score was produced with algorithm . The algorithm that most frequently produces high CE trees is  (with  very close behind) (see figs. <ref> and <ref> and table <ref>).

Let us now look at table <ref> and compare the results obtained for all algorithms with those of the real railway system and with the MST of the sites.
Algorithms ,  and  all consistently produced trees shorter than the railway, while all the trees produced by  were longer than the railway. The MST is significantly shorter than all the trees obtained by all algorithms.
All trees produced by all algorithms were more efficient than the MST. All algorithms produced trees with an on average larger efficiency (smaller $\nu$ value) than that of the railway. However, only for algorithm  were all trees produced more efficient than the railway.
Algorithms ,  and  all on average produced higher value CE trees than that of the railway or the MST. On the other hand, all trees produced by algorithm  showed CE values lower than that for the railway or MST.

It is important to note that a large CE value does not exactly equate to a better network. After all, the MST has a CE value 3% larger than the real railway tree, but it has a 31% larger $\nu$ value than the railway, making it much less efficient. CE was included in this work because it is a commonly used parameter in the literature, and in certain cases can help create a middle ground between construction cost and transport efficiency.

Overall, we can see that algorithm  maximizes transport efficiency (producing small $\nu$ trees), while algorithm  minimizes construction costs (producing small $L$ trees). However,  generates exceedingly long trees (its average tree presents a 21% difference to ) and is particularly inefficient computationally as mentioned previously, while 's average tree shows only a 4% difference in $\nu$ value compared to . Algorithms  and  show a middle ground between the two, being able to produce trees as efficient as  and slightly longer than  (their shortest trees show a less than 6% difference to 's shortest tree, and their average trees show a less than 2% difference to 's trees).

While all algorithms utilized were successful in producing more efficient trees on average than that of Portugal's mainland's railway system, we would overall recommend algorithm  when trying to produce an efficient (and somewhat short) tree, as it was the algorithm that showed, on average, higher CE values (better cost-efficiency) and lower $\nu$ values (higher efficiency), and also presented a smaller standard deviation with respect to $\nu$, thus more consistently produced efficient trees. This was similarly the same result obtained in <cit.>.

CHAPTER: CONCLUDING REMARKS

§ ACHIEVEMENTS

In this thesis, we aimed to test the adaptive Hagen-Poiseuille flow model (described in section <ref>) as a way to both mimic Physarum polycephalum's network adaptation dynamics and to be able to produce efficient networks.

We started by making a slight change to the model's implementation that allowed us to control the volume of fluid in each simulation, which is a key parameter that must be regulated to permit comparing between different runs of the simulation.

To determine whether or not the model studied produced realistic results and was worth studying, we tested certain key physical characteristics of the trees produced by it. We showed that the model tested was able to produce results that were consistent with physical observations, namely the tendency of the simulations to always choose the shortest path, especially as the continuous limit was approached. Additionally.
Murray's law (which is experimentally observed in several biological networks) was verified for bifurcation nodes in steady state trees, but it was not verified dynamically for the time of the simulation leading up to steady state.
By studying the results of the model applied to simple one-channel and two-channel networks, we observed that the conductivities and flows of said network's channels behaved like other biological contractile veins.

The model was also tested with regards to how well it could mimic Physarum polycephalum's network topology, using both fixed sources and sinks and non-static sources and sinks. Both cases were compared to Physarum's radial growth pattern, and both cases presented a key characteristic present in Physarum networks: showing fewer, thicker veins closer to the center of the network that branch out into many thinner veins as one gets closer to the periphery. However, the non-static sources and sinks steady states showed more accurate resemblance, as they also displayed other characteristics: connections between the sites on the periphery of the network (even though they were scarce compared to real Physarum networks) and thus the presence of loops in the network.
As such, it was determined that the model represents quite accurately Physarum-like network adaptation patterns, that is, it produces network shapes that resemble those of "older" (Physarum grows radially outward, so closer to the center of the network), more optimized parts of Physarum's networks. However, it is not successful at accurately describing the peripheral, foraging, "most recent" parts of Physarum's networks, that are highly connected by numerous thin veins and connect them using numerous loops.

The non-static sources and sinks configuration were also tested to determine whether or not they showed two key physical characteristics of Physarum: shuttle streaming and peristalsis. While this model does not aim to describe the chemical or biological dynamics of Physarum, it describes the basic hydrodynamics involved in the movement of cytoplasm inside its body. As this model has at its core the adaptation of vein walls of the network, it describes the contraction of the walls that is thought to contribute to shuttle streaming in Physarum.
It was shown that the simulations produced by the model presented peristalsis (oscillatory contraction and relaxation of the vein walls) on all edges tested.
It was also determined that the model did show shuttle streaming. This shuttle streaming was present in peripheral parts of the network (channels connecting sinks at the edge) and in parts of the network that were part of a loop, before a steady state was reached. As such, it is possible that the presence of loops in Physarum's trees may also contribute to the shuttle streaming phenomenon and the efficient transport of nutrients and resources throughout its cell.
Peristalsis and shuttle streaming were not observed for static sites configurations, and were only seen when using asynchronous sites. We can conclude that shuttle streaming and peristalsis are most likely consequences of the asynchronous resource consumption that occurs in Physarum, as mentioned previously.

Then, this thesis intended to determine if this model could be used to produce efficient networks, like Physarum has been witnessed to do; by efficient, one means networks with short distances between relevant sites. This was contrasted with determining whether or not the model was able to create short networks, like Steiner minimum trees.
Transport efficiency and network cost are two defining characteristics of a network that are heavily studied and weighed before a transport network is constructed. This biologically-inspired algorithm could be advantageous in constructing cost-efficient networks if proven to be successful.
The implementation of the model involves a stochastic algorithm, as the final shape of a network produced by this model is heavily influenced by the starting conditions of conductivities, which are set to be random; like Physarum, this model can produce many different types of networks, depending on its environment. The choice of sources and sinks was stochastic as well and followed one of four different tested algorithms. The results were tested on two different types of networks: first, connecting $n=\{3,4,5\}$ sites in a regular polygon configuration; then, connecting $n=25$ sites in a configuration that is meant to mirror mainland Portugal, thus creating a realistic transport network.
All algorithms collectively were shown to be able to produce short trees, efficient trees, and cost-effective trees, and all produced trees that were more efficient than the studied Portuguese railway network. Algorithm  produced overall long, quite efficient trees (about 20% longer than  for the Portuguese configuration), and was the most computationally demanding one. Algorithm  produced mainly short, less efficient trees (yet showing only about a 4% difference in efficiency to ). The recommended algorithm to use to produce cost-efficient trees would be , as it presented the overall highest efficiency values (behind ) with reduced length values and adequate computation time.

§ FUTURE WORK

The model studied in this work does not fully mirror Physarum's network characteristics; namely, it doesn't mimic the peripheral, "more recent" parts of the trees that often present numerous thin veins connected by many loops. However, the model produces networks that share similarities with the "older", more optimized parts of Physarum's network. Thus, the model should be extended to include a description of the growth and foraging mechanisms of slime mold that can accurately reproduce network patterns.

While shuttle streaming was observed using this model, it wasn't observed at steady state or in all edges tested. The model used does not explicitly incorporate peristaltic mechanisms <cit.> or other mechanisms that are theorized to be responsible for shuttle streaming like the transport of signaling molecules released in regions that come in contact with stimuli <cit.>. Including such mechanisms in the model could result in the active presence of shuttle streaming and more accurate network patterns.

The general form of the adaptation equation (eq. (<ref>)) allows for any choice of $g$ function. In this work, only the $g$ function described in eq. (<ref>) was used, as it minimizes the dissipated power. It would be interesting to study the results of using a different $g$ function, particularly one that saturates for large flux values to mimic what is observed in real biological networks (like blood vessels <cit.>).

The study of the efficiency of the networks created using this model could be extended by using different $g$ functions (like $g = |Q_{ij}|^\gamma$ for different $\gamma$ values, like as seen in <cit.>) to create them, as this could lead to different topologies and efficiencies. This study could also be significantly improved by optimizing the model's computational implementation to reduce run times; as this is a stochastic algorithm that relies on running the simulations numerous times, optimizing the run time of the algorithm (by, for example, parallelizing the code) is essential.

There's still much work to be done in modeling all of Physarum polycephalum's complex behavior. After all, we still lack a full understanding of the mechanisms behind its cytoplasmic movement, among other aspects. Because of this, we have yet to find a consolidated model that can mimic Physarum's network dynamics and intelligent behavior. This is an ambitious open-ended problem that could have exciting and beneficial consequences for a large range of fields.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

howpublishedSlime Mold - Georgios Cherouvim, <https://ch3.gr/slime-mold/>
year2018, accessed in May of 2023

American Physical Society

journaltitlePhys. Rev. Lett.
titleDamage and Fluctuations Induce Loops in Optimal Transport Networks

Society for Industrial
Applied Mathematics

titleMathematical Problems in the Biological Sciences

journaltitlePhysica A: Statistical Mechanics and its Applications
titlePhysarum solver: A biologically inspired method of road-network navigation

journaltitlePhysical review letters
titleAdaptation and Optimization of Biological Transport Networks

journaltitlePhysical review letters
titleStructure, Scaling, and Phase Transition in the Optimal Transport Network

journaltitlePhysical Review E
titleDirectional sensing and streaming in Dictyostelium aggregation

titleA Method Inspired by Physarum for Solving the Steiner Problem

journaltitleRoyal Society Open Science
titleRemarkable problem-solving ability of unicellular amoeboid organism and its mechanism

Berlin, Heidelberg

Springer Berlin Heidelberg

booktitleProbability on Discrete Structures
titleModels of First-Passage Percolation

titleMaze-solving by an amoeboid organism

titleRules for Biologically Inspired Adaptive Network Design

titleAdvances in Physarum Machines: Sensing and Computing with Slime Mould
family=Le Bernemi,

howpublishedWikipedia, <https://commons.wikimedia.org/wiki/File:Blob_(Physarum_polycephalum).jpg>
titleBlob (Physarum polycephalum)
year2020, accessed in May of 2023

howpublishedWikipedia, <https://en.wikipedia.org/wiki/File:Physarum_polycephalum_plasmodium.jpg> (originally published on flickr)
titlePhysarum polycephalum plasmodium
year2011, accessed in May of 2023

journaltitleJournal of Physics D: Applied Physics
titlePhysarum polycephalum—a new take on a classic model system
family=Tim Tim (VD fr),

howpublishedWikipedia, <https://commons.wikimedia.org/wiki/File:Physarum_polycephalum_exploration_2.jpg>
titlePhysarum polycephalum exploration 2
year2021, accessed in May of 2023

howpublishedThe Slime Mould Collective, <https://slimoco.ning.com/photo/physarum-polycephalum-on-agar?context=user>
titlePhysarum Polycephalum on Agar
year2015, accessed in May of 2023

journaltitleProceedings of the National Academy of Sciences
titleRandom network peristalsis in Physarum polycephalum organizes fluid flows across an individual

journaltitleProceedings of the National Academy of Sciences
titleMechanism of signal propagation in Physarum polycephalum

journaltitlePhysical review letters
titleAmoebae anticipate periodic events

journaltitleProceedings of the National Academy of Sciences
titleSlime mold uses an externalized spatial "memory" to navigate in complex environments

journaltitleBehavioral Ecology
titleAmoeboid organism uses extracellular secretions to make smart foraging decisions

titlePhysarum-inspired Network Optimization: A Review

journaltitlePhysics of Life Reviews
titleDoes being multi-headed make you better at solving problems? A survey of Physarum-based models and computations

Massachusetts Institute of Technology Press (MIT Press)

journaltitleArtificial Life
titleCharacteristics of pattern formation and evolution in approximations of physarum transport networks

journaltitleJournal of theoretical biology
titleA mathematical model for adaptive transport network in path finding by true slime mold

journaltitlePhysica D: Nonlinear Phenomena
titleAdaptive Hagen–Poiseuille flows on graphs

titleFormation and Optimisation of Vein Networks in Physarum

journaltitleJournal of Theoretical Biology
titleFlow of a viscous fluid through an elastic tube with applications to blood flow

journaltitleProceedings of the National Academy of Sciences of the United States of America
titleThe Physiological Principle of Minimum Work: I. The Vascular System and the Cost of Blood Volume.

titleA Method Inspired by Physarum for Solving the Steiner Problem

journaltitleNetworking, IEEE/ACM Transactions on
titleOn the Solution of the Steiner Tree NP-Hard Problem via Physarum BioNetwork

journaltitleRoyal Society Open Science
titleRemarkable problem-solving ability of unicellular amoeboid organism and its mechanism

journaltitleJournal of Experimental Biology
titleOscillatory Contraction Activity in Physarum

International Press of Boston

journaltitleCommunications in Mathematical Sciences
titleODE- and PDE-based modeling of biological transportation networks

journaltitleProceedings of the National Academy of Sciences
titleThe Physiological Principle of Minimum Work

journaltitleThe Journal of general physiology
titleOn connecting large vessels to small. The meaning of Murray's law

journaltitleJournal of Physics D: Applied Physics
titleExperimental models for Murray’s law

howpublishedPerlebär - Martins Seiten aus der Prignitz, <https://prlbr.de/2019/euclidean-steiner-trees-in-regular-polygons/>
titleEuclidean Steiner trees in regular polygons
year2019, accessed in May of 2023

Boston, MA

Springer US

titleReducibility among Combinatorial Problems

journaltitleBiophysical Chemistry
titleSmart network solutions in an amoeboid organism

booktitle2012 IEEE Wireless Communications and Networking Conference (WCNC)
titleA physarum-inspired algorithm for minimal exposure problem in wireless sensor networks

journaltitleBMC Systems Biology
titleA physarum-inspired prize-collecting steiner tree approach to identify subnetworks for drug repositioning

booktitle2016 IEEE Congress on Evolutionary Computation (CEC)
titleFast algorithms inspired by Physarum polycephalum for node weighted steiner tree problem with multiple terminals

journaltitleBMC Bioinformatics
titleThe node-weighted Steiner tree approach to identify elements of cancer-related signaling pathways

titleBiological flow networks: The absolute basics
year2019, accessed in May of 2023

journaltitleInt J Bifurcat Chaos
titleSlime Mould Inspired Generalised Voronoi Diagrams with Repulsive Fields

The MIT Press

booktitleIntroduction to Algorithms, Fourth Edition
titleBreadth-first search

CHAPTER: STEINER TREES ADDITIONAL IMAGES

§ TRIANGULAR CONFIGURATION

§ SQUARE CONFIGURATION

§ PENTAGONAL CONFIGURATION

§ MAINLAND PORTUGAL'S CONFIGURATION