A generative modeling approach to connectivity—Electrical conduction in vascular networks

Journal of Theoretical Biology 399 (2016) 1–12

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A generative modeling approach to connectivity—Electrical conduction in vascular networks Bjørn Olav Hald Department of Biomedical Sciences, University of Copenhagen, Copenhagen, Denmark

H I G H L I G H T S

G R A P H I C A L

A B S T R A C T


Strategy for dynamical models: Generate connectivity using patternalgorithms.
Software generates electrical models of arbitrary vascular networks from cell level.
Series of electrical conduction models generated in various vascular networks.
Vascular morphology and connectivity delimit conduction capacity within networks.

art ic l e i nf o

a b s t r a c t

Article history: Received 25 November 2015 Received in revised form 7 March 2016 Accepted 18 March 2016 Available online 31 March 2016

The physiology of biological structures is inherently dynamic and emerges from the interaction and assembly of large collections of small entities. The extent of coupled entities complicates modeling and increases computational load. Here, microvascular networks are used to present a novel generative approach to connectivity based on the observation that biological organization is hierarchical and composed of a limited set of building blocks, i.e. a vascular network consists of blood vessels which in turn are composed by one or more cell types. Fast electrical communication is crucial to synchronize vessel tone across the vast distances within a network. We hypothesize that electrical conduction capacity is delimited by the size of vascular structures and connectivity of the network. Generation and simulation of series of dynamical models of electrical spread within vascular networks of different size and composition showed that (1) Conduction is enhanced in models harboring long and thin endothelial cells that couple preferentially along the longitudinal axis. (2) Conduction across a branch point depends on endothelial connectivity between branches. (3) Low connectivity sub-networks are more sensitive to electrical perturbations. In summary, the capacity for electrical signaling in microvascular networks is strongly shaped by the morphology and connectivity of vascular (particularly endothelial) cells. While the presented software can be used by itself or as a starting point for more sophisticated models of vascular dynamics, the generative approach can be applied to other biological systems, e.g. nervous tissue, the lymphatics, or the biliary system. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Connectivity generation Modeling Vascular networks Electrical conduction

1. Introduction Abbreviation: EC, endothelial cell; GJ, gap junction; MEGJ, myoendothelial gap junction; ODEs, ordinary differential equations; SMC, smooth muscle cell; Vm , membrane potential; Vm;rest , resting membrane potential E-mail address: [email protected] http://dx.doi.org/10.1016/j.jtbi.2016.03.032 0022-5193/& 2016 Elsevier Ltd. All rights reserved.

The integrative behavior of a biological system depends on the dynamic interactions among its constituent parts. Whereas reductionism advocates study of constituent parts, biological

2

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

Fig. 1. Basis of generative modeling approach to connectivity. (A) Vascular beds differ structurally. Structural differences are obvious from network topology but can also exist at both vessel and cell levels. The basic idea is that larger structures are generated from smaller structures via “pattern generators” or “tabulated data” (see text). Usually, tabulated data is more frequent at macroscopic levels and pattern generation occurs at microscopic levels. (B) The “lego block principle” uses object-oriented programming to describe each entity (only low level entities are depicted) of the biological system as a type that can hold a unique index and entity specific parameters. Sets of ODEs describing (1) internal dynamics of each entity, and (2) interactions among the entities are collected. Using a pattern generator algorithm and/or tabulated data, the calculated number of entity-objects can be copied from each entity-type and the relevant sets of connectivity produced, one set for each type of interactions. Objects can be arbitrarily organized but connectivity generation must be based on biological observations.

function usually arises from multicellular interactions. Indeed, many biological systems display salient structural features that influence global behavior, e.g. the network structure of different vascular beds (Cassot et al., 2006; Hirsch, 2012; Segal, 2005), see Fig. 1A. The isolated impact of structure can be studied using mathematical modeling. Modeling of structural changes in spatially extensive systems, however, is laborious as changes in connectivity necessitate model reorganization. Here, a general strategy to ease connectivity generation across multiple spatial scales in dynamical models is presented using electrical signaling in the vascular wall as an example. Microvascular networks represent highly dynamical and extensive spatial systems that rapidly redistribute blood flow to underperfused areas within a tissue (Hill, 2012; Segal, 2005). Networks display distinct structural differences between vascular beds, e.g. skeletal and cerebral circulations (Segal, 2005; Cassot et al., 2006), yet the structural impact on vascular tone regulation is relatively unknown (Hirsch, 2012). Conversely, constituent endothelial cells (ECs) (Nilius and Droogmans, 2001; Mehta and Malik, 2006) and smooth muscle cells (SMCs) (Longden et al., 2016; Yamamoto et al., 2001; Cole and Welsh, 2011; Hill et al., 2006) are well-studied as they underlie a host of regulatory mechanisms controlling local vascular tone (Segal, 2005; Somlyo et al., 1999; Hill et al., 2006). Electrical communication between vascular cells is critical to integration of the local vasomotor responses across a network (Tran et al., 2012; Hill, 2012; Segal, 2005). A local electrical signal rapidly disperses via gap junctions (GJs) between vascular cells. Electromechanical coupling in SMCs converts the signal into a conducted vasomotor response that coordinates tone across long distances within the network (Emerson et al., 2002; Tran et al., 2012; Xia and Duling, 1995; Segal et al., 1989). Theoretical studies have shown that not only electrical properties of ECs and SMCs but also cell orientation and number of SMC layers (i.e. capacitive load) influence conduction (Diep et al., 2005; Kapela et al., 2010; Hald et al., 2015). We therefore hypothesize that connectivity and morphology of vascular structures in

general delimit electrical communication within a vascular network. Consequently, the aims of the present paper were threefold: (a) To present a novel strategy to generate connectivity in computational models of dynamical systems. (b) To apply the strategy in an open source software that generates dynamical models of arbitrary vascular networks from cell level to the network level (in the following, the software is called ‘the vascular model generator’). Each level may be furnished with particular dynamics and/or connectivity. (c) To use the vascular model generator to study the impact of network structure and topology on electrical spread (see below). Existing software platforms mostly focus on dynamics and changes herein, notably biochemical reaction-networks, intracellular cell-signaling or regulatory systems (see http://www.sys tems-biology.org/software/simulation/). Except from large projects such as ‘TheVirtualHeart.org’, few platforms handle spatial changes (Stoma et al., 2011; Cornelis et al., 2012) and none addresses spatially extensive systems. Here, a “lego block principle” is used to model the basic entities of a system (the EC, the SMC, and the vessel) and their interactions. Subsequently, connectivity between entities is generated from either algorithms that approximate the organization of entities or from tabulated morphometric data. The presented vascular model generator can therefore be used in conjunction with recent methods for detailed 3D reconstructions of vascular network topology based on integrative 3D imaging from disparate sources (Rieger and Welter, 2015; Kim et al., 2012; Blinder et al., 2013; Hirsch, 2012), ex vivo vascular castings (Cassot et al., 2006; Hirsch, 2012), or synthetic reconstructions (Schneider et al., 2014). Such reconstructions are commonly used to analyze network topology (Blinder et al., 2013; Cassot et al., 2006; Hirsch, 2012) or to model e.g. hemodynamics (Kim et al., 2012; Rieger and Welter, 2015; Reichold et al., 2009) or vascular remodeling (Pries and Secomb, 2014; Jacobsen et al., 2003). Due to the separation of dynamics and connectivity, the vascular model generator can also be applied to model dynamics of other phenomena in vascular networks, e.g. nephron synchronization or vascular remodeling

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

3

(Pries and Secomb, 2014; Marsh et al., 2013). Here, simple dynamical models of electrical behavior within vascular cells were used to generate structurally different models of electrical conduction along the vascular wall. In keeping with the our aims, we specifically asked: to what extent is electrical spread along the vascular wall affected by (a) connectivity and conductance of GJs, (b) the size of vascular cells, and (c) the connectivity of vessel segments within the network?

2. Materials and methods The present study extends our previously presented strategy (Hald et al., 2013) that, instead of building a dedicated program for a specific single model, emphasizes development of a model generator program that produces computational models to be simulated based on separation of structural from kinetic input. Basically, living systems are generally organized in hierarchies via more or less regular assemblies of a limited set of smaller entities, e.g. a vascular network consists of blood vessels that in turn are composed of ECs and SMCs, see Fig. 1A. All biological entities (structures or compartments) that constitute the model are identified and each entity is modeled as a type using object-oriented programming (OOP) (Hald et al., 2013). The sets of ordinary differential equations (ODEs) that describe dynamics (1) within each entity and (2) between entities are defined separately (Fig. 1B). Once defined, the entities can be connected in arbitrary ways to build a biological structure (lego block principle); this constitutes the basis of model generation. In general, connectivity generation proceeds as follows: (1) Copies of each type of entity of the system are stored as objects each with a unique ID and possibly unique model parameters. Each set of same-type-objects is associated with the corresponding set of ODEs defining internal kinetics for the type. (2) The organization of objects can almost always be approximated by some connectivity generating algorithm (pattern generator) or be defined by tabulated (external) data. Thus, multiple sets of connectivity, with each set describing all couplings between two types of entities, are produced. The spatial component in solving a dynamical system can be dissolved by only addressing dynamic exchanges between objects. Thus, each set of connectivity is associated with a relevant set of ODEs defining dynamics of interactions. Use of OOP and separation of dynamics from connectivity generation therefore allows for construction of more or less complicated models of e.g. vascular networks based on either synthetic or real morphometric data. 2.1. The vascular model generator A vascular network is conveniently represented by a directed graph of edges (vessel segments) and nodes (internal branch points and network boundaries) (Reichold et al., 2009). Each vessel is composed of an EC layer circumscribed by one or more SMC layers. SMCs are oriented in a circumferential direction, perpendicular to ECs (see “Vascular network” in Fig. 2). A network graph is constructed from the set of 3D node coordinates delineating vessel segments. Each vessel is assigned a resting diameter, length, number of cell layers, and the average size of ECs and SMCs within the layer (if not assigned, default values are used). Thus, the network graph is constructed from tabulated data (that also can be generated, e.g. from synthetic network reconstructions). The network information is used to calculate the number of all modeled objects based on network, vessel and cell sizes. A set of patterngenerator algorithms (see below) establishes cell–cell connectivity. The information is stored in a hierarchy that follows the hierarchy of modeled structures (gray elements of “Data hierarchy” in Fig. 2). The vascular model generator program is written in Python 2.7.9

Fig. 2. Structural and connectivity information obtained for a specific vascular network is converted into a data hierarchy that organizes the system of variables and ODEs. (A) Any “real” vascular network can be converted into a network graph either (1) by use of kinks to straighten tortuous vessel segments or (2) by directly specifying vessel length (overriding default calculation of vessel length between nodes). Each edge/vessel (see insert), minimally described through node coordinates (2D or 3D), is passed to the model generator that calculates connectivity and numbers of cells within cell layers according to a preselected pattern (see inset bottom). Subsequently, connectivity between vessels are calculated. (B) Connectivity of all biological structures (gray boxes) in a vascular network is stored in a data hierarchy as well as structure specific parameters. Thus, the toplevel network object stores NV vessel objects. Each vessel object stores 1 EC layer, and NS SMC layer objects. Each cell layer stores a number of cells (NEC or NSMC ) as well as the cell–cell coupling vector(s) (white boxes) associated with gap junctional currents. (C) The data hierarchy organizes all parameters and indices used to simulate the evolution of variables stored in vector, y, i.e. the organization of y follows straightforwardly from the data hierarchy.

(files available in the Supplemental information) allowing for rapid prototyping, e.g. alteration of pattern-generator algorithms. As described in Hald et al. (2013), the vascular model generator produces (1) a “connectivity data-file” and (2) a “C þ þ model-tobe-simulated” (for fast execution of simulations). Briefly, the C þ þ program uses the data-file to establish an extended data hierarchy containing all data that can be set or calculated before the ODEs are solved. This hierarchy unfolds as a natural traversal of every cell object within the network. As each cell accesses its own parameters and connectivity information, all relevant variables (see “Variable vector” in Fig. 2) and parameter values are available to calculate y_ ¼ f ðyÞ (the right-hand side of all ODEs) and/or Jacobian information without storing of matrices. Using a finitedifference method, the set of ODEs (initial value problem) is solved by the GMRES solver from CVODE-2.5.0 (Sundials) with a relative

4

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

tolerance of 10  4 and absolute tolerance of 10  9. The C þ þ program was compiled using GNU gcc-4.5.0. 2.2. Electrophysiological dynamics Due to the focus on structural modeling, only basic electrical dynamics within and between ECs and SMCs are considered in the model. Thus, membrane potential ðVm Þ is determined by the cell capacitance and the sum of currents (a) across the cell membrane (internal dynamics) and (b) through GJs between neighboring cells (interactions) (Diep et al., 2005). All currents are assumed to be proportional with the size of relevant cell surface area, directly linking morphology to electrophysiology. Experimental measures of cell size and cellular dimensions within arterioles display large standard deviations and average sizes change between vascular beds and/or species as well as during development (Haas and Duling, 1997; Sandow et al., 2003a,b). Here, cell sizes were simply approximated through rectangular dimensions, length (l), width (w), and (radial) height (h), see Fig. 4A, and based on literature averages: EC: 100  7  0:75 μm3 (Haas and Duling, 1997; Sandow et al., 2003a, 2003b) and SMC: 75  7  5 μm3 (Haas and Duling, 1997; Gattone et al., 1986). Current across the membrane can be measured in physiological solutions using the whole-cell patch-clamp technique and may thus be described through a simple non-linear resistor, i.e. an IVcurve, represented via a function Im ðVm Þ. As current density is normally measured, we assume a specific membrane capacitance of c ¼ 1 μF=cm2 and that current therefore scales with the cell surface area, here, approximated by Am ¼ 2  ðhw þ hl þwlÞ. The current density functions for an EC and a SMC, respectively, were fitted from Schnitzler et al. (2000) and Jackson et al. (1997) (see Fig. 3): 2 IEC  V2m þ 1:00  10  4  V3m m ðVm Þ ¼ 19:59 þ 0:99  Vm þ 1:66  10 2 SMC Im ðVm Þ ¼ 35:58 þ 1:68  Vm þ 2:48  10  V2m þ 1:29  10  4  V3m

Note that vascular cells are voltage-gated, but non-excitable (Hill, 2012; Wolfle et al., 2011; Hald et al., 2014). Contrarily, currents through GJs are assumed to be ohmic under physiological conditions ðΔVm 7 20 mVÞ (Vogel and Weingart, 2002). Homocellular couplings (homo) are located along cell edges and myoendothelial couplings (mec) are found between cell layers (see Table 1), i.e. GJ conductance (gGJ) is assumed to scale with border area approximated by Ahomo ¼ 2h  ðw þ lÞ and Amec ¼ w  l, A ⋆ respectively: g GJ ¼ AGJ ⋆  g GJ , where GJ represent homocellular or

myoendothelial gap junctions (MEGJs) and ⋆ denotes standard value. Standard values are given in Table 1 (Yamamoto et al., 2001; Lidington et al., 2000). With superscript referring to individual cells, the membrane potential ðVim Þ for cell i therefore evolves by: C im

dVim ¼ C im  I m ðVim Þ Iistim dt þ

N X Ahomo g ⋆

homo

j

þ



Ahomo

ðinternal cell dynamicsÞ

Vim  Vjm

M   X Amec g ⋆ mec Vim  Vkm A mec k



ðhomocellular interactionsÞ

ðmyoendothelial interactionsÞ ð1Þ

where j designates each of the N cells coupled via homocellular GJs, and k designates each of the M cells coupled via MEGJs. Istim denotes current stimulation input (zero by default). As a large vascular network inherently generates a large, sparse coupling matrix, a matrix-free solver to the system V_ m ¼ f ðVm Þ is employed (Hald et al., 2013). This, in turn, entails a proper handling of cellular connections. GJ currents represent the only interactions between cells, i.e. currents are calculated using the relevant sets of connectivity (see Figs. 1 and 2 – for more detail, see Section S3 of the Supplemental Information).

Table 1 Reference GJ conductances and calculated surface areas. As gcell depends on the specific GJ conductance (conductance per surface area) and the relevant surface area of the cell, it is calculated as gcell ¼ A  g⋆ =A⋆ , where ⋆ denotes the standard cell value. The conductance for each individual GJ plaque between two cells (gGJ) depends also on the number of couplings per cell (N), i.e. gGJ ¼ g cell =N. Cell type Otherwise

EC SMC

g⋆ homo;EC

g⋆ homo;SMC

g⋆ mec

ð3 MΩÞ  1 2 A⋆ m ðμm Þ

ð90 MΩÞ  1 2 A⋆ homo ðμm Þ

2 A⋆ mec ðμm Þ

1560 1870

160 820

700 525

ð900 MΩÞ  1

GJ

Fig. 3. Vm within a cell depends on membrane currents and gap junction currents. (A) The current densities of every EC (gray curve) and SMC (black curve) used in simulations. Vm;rest ð
Þ is  40 mV for both cell types. If Vm o Vm;rest , outward current drives Vm back toward Vm;rest and vice versa (see arrows in A). (B) Total coupling conductance per cell (gcell) is treated as an extensive property. Couplings between cells are facilitated by GJs (black arrows through pores) embedded in the plasma membrane of a cell (gray shade). The relevant membrane area for GJs is denoted ‘A’ and the number of couplings per cell is denoted ‘N’. As cell b is larger than cell a, b harbors more GJs and g cell;b 4 gcell;a . Cells b and c have equal areas, i.e. gcell;b ¼ g cell;c , but as c has more couplings, coupling strength per coupling in c is lower as each coupling must share gcell;c .

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

5

Fig. 4. Auto-generation of network connectivity is based on the physical dimensions of vessel segments and their constituent cells. Upper panel: Vessel connectivity. Cell length and width are denoted by lx and wx with x representing either SMC or EC. (A) To calculate the approximate number of cells in a vessel segment, the cell height (radial dimension) is neglected and each cell is assumed to be of rectangular shape whereas the vessel segment is assumed to be cylindrical with length, lv, and a resting diameter, dv. (B) ECs or SMCs are then filled into approximately π  dv  lv rectangular grids (due to the discrete nature of vascular dimensions, exact matches are obviously not possible). General algorithms establish cell indices as well as homocellular and myoendothelial connectivity. These algorithms only require the number of cells in axial and circumferential directions and a chosen number of couplings (illustrated by arrows). (C) By connecting cells at the longitudinal grid borders, a tubular vessel connectivity is established. Lower panel: Branch connectivity. (D) Due to morphology, branch point connectivity in the EC and SMC layers is different. In both cases, the process starts by aligning cells within the last and first circumference, respectively, of connecting branches. (E) and (F) EC layer connectivity: cells of each branch are partitioned in sequential groups, the size of each corresponds to the weighted sum of cells in the circumference of other branches (see color coding). Each same-colored partition of each branch is coupled such that each individual cell attains its designated number of couplings per cell. SMC layer connectivity: cells of each branch are distributed in a one-to-one fashion to cells within other branches: each SMC couples to (at least) one SMC from each other branch. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

2.3. Connectivity generating algorithms Pattern-generating algorithms calculate cell connectivity within 2D or 3D networks based on a few simple ideas (Fig. 1A): First, a set of simple patterns establishes how cells are distributed and coupled within each vessel segment. The vessel segment is assumed to be cylindrical and comprised of rectangular ECs and SMCs. Second, a set of simple patterns establishes how cells couple across branch points, i.e. between vessel segments. Using the two algorithms, connectivity of cells within any vascular network described via a graph can be calculated. To utilize available network reconstruction data, a graph is described via tabulated data. The network graph represents connectivity between vessel segments. These ideas, from a single vessel segment to a full network, are unfolded in the following (refer to Fig. 4 throughout).

2.3.1. Generation of cell connectivity in a single vessel segment A cylindrical vessel segment has length (lv) and resting outer wall diameter (dv).

(a) The average length, width, and height ðl  w  hÞ of ECs or SMCs within a vessel segment are used to calculate the approximate number of cells in the circumferential (ncirc) and longitudinal (naxis) directions in each cell layer (Fig. 4A): d  π  ncirc;SMC  lSMC  ncirc;EC  wEC

ð2Þ

l  naxis;SMC  wSMC  naxis;EC  lEC

ð3Þ

Strict equality between the two sets of equations is dismissed because of the elasticity of wall-material and the plasticity of vascular cells that allows vessel segments to dilate or contract substantially without concomitant increase in the number of cells, ncells ¼ naxis  ncirc for each cell layer. Each cell is assigned a unique integer index. The rectangular shape eases automatic assignment of connectivity but cells are obviously not brick shaped. Importantly, the shape or size of a cell is therefore not directly coupled to connectivity. This allows for morphological heterogeneity between cells (e.g. membrane area or cell volume). (b) Homocellular cell–cell couplings are calculated by choosing a particular coupling pattern (e.g. a 4 neighbor pattern can be

6

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

calculated via nj ¼ ni 7 k; k A ½1; ncirc , where i and j represent the current and neighboring cells, respectively). Here, we used either 4, 6 or 8 couplings for each cell (see red arrows in Fig. 4B). For all types of couplings, a connection is only stored once, i.e. if cells i and j are coupled, only a single connection ½i; j is stored. (c) Myoendothelial connectivity, i.e. EC:SMC couplings, is slightly more complex given the perpendicular orientation of ECs and SMCs. Usually, the EC and SMC “grids” do not superimpose ncirc;EC n evenly (i.e. ncirc;SMC 2 = N and/or naxis;SMC 2 = N). A random generator is axis;EC then used to determine which cells have more (or less) connections to avoid introduction of systematic errors. Furthermore, heterogeneity both in the number and location of EC: SMC couplings can be introduced by choosing a random subset of connections (again, using a random number generator) to be used in the model. The procedures are encapsulated in a set of pattern-generating algorithms that together determine connectivity within each cell layer of a single vessel segment. Connectivity is changed simply by altering one or more of these algorithms.

2.3.2. Generation of vessel connectivity within the network Vessel segments are successively incorporated into the network. The first cell of each new vessel is given an index offset equal to the current sum of cells in the current network. Thus, each cell has a unique index, allowing for connectivity across branch points to be calculated in the following manner (Fig. 4D to F): (a) As the network graph is directed, the low and high index ends of each cell layer of each vessel segment are known. (b) Each branch point/internal node of the graph is traversed (outer loop). (c) In each branch point (inner loop), the circumference of cells from each cell layer of “incoming” and “outgoing” branches are collected (Fig. 4D–F). For each layer, these “branch cells” can be connected in various ways (see Section S4 in the Supplemental Information). Here, by default, Branch SMCs are distributed in an “all-to-all” fashion: due to the circumferential orientation, every SMC from each branch couples to a single SMC from all other branches. Branch ECs are divided according to the number of branches: each partition of each branch holds a weighted number of ECs in proportion to the relative number of nEC circ in other branches. For example, if one small and two large branches connect, then most ECs from a large branch will couple to ECs from the other large branch and only few ECs couple to ECs from the small branch. The number of EC:EC couplings per cell at the branch point equals EC:EC couplings per cell within vessel segments. (d) Couplings from each branch point are collected in sets of connectivity, one set for each cell type (EC or SMC) that connects. Following completion of connectivity generation, the sets of connectivity are written to disk and used by the C þ þ simulation program. As previously described (Hald et al., 2013), connectivity and model parameters are assembled in a data hierarchy (Fig. 2) that is traversed upon calculation of V_ m ¼ f ðVm Þ. The unique index of each cell i in the network corresponds to the location of Vim in vector Vm . During a simulation, Vm is stored at each time step. Hence, each Vim variable follows the same indexing as used in the connectivity data-file that organizes the data hierarchy. This simplifies data extraction for visualizations.

2.4. Simulations: electrical spread in vascular networks Any cell can be electrically stimulated either by (1) injecting current ðIstim a 0Þ or by (2) applying a voltage clamp (setting V_ m ¼ 0). Stimulation elicits a current between coupled cells (see Eq. (1)), i.e. an electrical conduction process, that spreads along the vascular wall. In simulations, electrical spread along the vascular wall was initiated by injection of hyperpolarizing current (i.e. Iistim o 0 pA) or by a hyperpolarizing voltage clamp (setting Vim o  40 mV) in one or more ECs for 150 ms. To illustrate how morphology, connectivity, and topology affect electrical communication within the vascular wall, either connectivity, cell dimensions, or the vessel topology was subject to systematic changes between simulations. In each model, every EC and every SMC, respectively, were identical and all kinetic parameters remained constant across simulations.

3. Results Structurally different models of electrical signaling were generated at three levels: within a single vessel, across a single branch point, and within vascular networks. 3.1. At the single vessel level Most modeling studies on electrical conduction use parameter values of GJ conductances that disregard assigned number of couplings per cell and/or cell sizes (Diep et al., 2005; Hald et al., 2012; Kapela et al., 2010; Wolfle et al., 2011). The GJ conductances used are based on a limited supply of rough experimental estimations (Lidington et al., 2000; Yamamoto et al., 2001; Rijen et al., 1997). Here, coupling conductance is treated as an extensive property and therefore scaled to the number of couplings and the relevant cell surface area (see Table 1). Number of couplings: The EC layer is electrically well-coupled and composed of longitudinally oriented cells. These properties are critical in long-range conduction of electrical perturbations (Diep et al., 2005; Hald et al., 2012; Tran et al., 2012; Wolfle et al., 2011). First, the dependence of conduction on the number and connectivity of EC:EC couplings was investigated – all other parameters remained constant between simulations (see Fig. 5A). If EC:EC conductance per (individual) GJ is a constant, an increase in the number of couplings also increases the overall coupling strength of a cell and, in turn, of the entire vessel segment. This is shown in Fig. 5B, where an increasing number of couplings causes less conduction decay as well as smaller local Vm -responses (due to lower input resistance). Conversely, if EC:EC conductance per cell is treated as an extensive property, the overall coupling strength becomes independent of the number of couplings (see Table 1). As shown in Fig. 5C (fully drawn curves), however, more couplings still cause a smaller conduction decay along the vessel wall. This is a consequence of additional couplings being recruited from diagonally located ECs (see Fig. 4B), i.e. more couplings give rise to a larger proportion of GJs being longitudinally oriented. Thus, a coupling pattern that favors longitudinal couplings produces less conduction decay compared to more circumferential or radially oriented coupling patterns. Finally, using a scaled EC:EC conductance of 13 μS per cell (dashed line in Fig. 5C), the electrical conduction decay is higher than experimentally observed (Xia and Duling, 1995; Segal et al., 1989; Emerson et al., 2002) (compare dashed curve to the gray shaded area, showing experimental decay constants of  1–2 mm using the short cable equation Hald et al., 2012). Therefore, specific GJ conductances (g ⋆ GJ per cell) were multiplied by 8 in all subsequent simulations to align conduction decay with experimental observations.

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

7

Fig. 5. Single vessel level: connectivity affects electrical communication. (A) Vessels differing only with respect to the number of couplings from an EC to neighboring ECs were constructed. A hyperpolarizing current injection of 900 pA into each ECs at the vessel end elicited an electrical spread that was assessed by measuring Vm in the SMC layer along the vessel axis. (B) GJ conductance was set to 1=ð3 MΩÞ per coupling. Each EC has either 4 (black), 6 (gray), or 8 (light gray) couplings, i.e. total conductance per cell increases with number of couplings. (C) Solid curves: GJ conductance was set to 8=ð3 MΩÞ per cell, i.e. conductance per EC:EC coupling decreases with number of couplings. Dashed curve: GJ conductance was set to 1=ð3 MΩÞ per cell. Gray shaded area: exponential decay range using the short cable equation, 1000 μm o λ o 2000 μm.

Due to the circumferential orientation of SMCs and their relatively low GJ conductance (Diep et al., 2005; Hald et al., 2012; Tran et al., 2012; Hill, 2012; Wolfle et al., 2011), physiologically unrealistic increases in the number of SMC couplings were required to affect longitudinal electrical spread (e.g. using the same protocol as in Fig. 5B, the Vm -difference between 4 and 8 SMC:SMC couplings was less than 0.05 mV at 2000 μm). Changing cell dimensions: In the set of simulations shown in Fig. 6, the effect of changing (1) EC length ð 725 μmÞ, (2) EC width ð 7 2 μmÞ, or (3) SMC width ð 7 2 μmÞ is addressed (all other parameters were held constant between simulations; the number of couplings per cell is set to six). Smaller cells cause a limited increase in membrane area and overall capacitance (for ECs o 8%, for SMCs o16%) as cell height (radial dimension) is kept constant. The surface area of each cell layer, however, is constant, i.e. overall MEGJ conductance is constant. In Fig. 6B, conduction was initiated by injection of  900 pA current into local ECs which produces local Vm -differences because total current injection is proportional to the number of circumferential ECs. This effect is mitigated using a voltage clamp ð  55 mVÞ as stimulus in local ECs (Fig. 6C). Particularly EC length and thickness affect electrical communication along the vascular wall. Longer ECs result in fewer membranes to cross in the longitudinal direction and more membrane area with potential GJ plaques. Hence, electrical spread is enhanced. Slimmer ECs result in more cells in the circumferential direction per vessel and an increase in overall membrane area for GJs. Thus, overall conductance increases, causing enhanced electrical spread. In addition, more circumferential cells give rise to more locally stimulated cells, i.e. more charge is available for longitudinal spread. Changes in SMC thickness also changes the number of SMCs (but not MEGJ conductance along the vessel wall). However, slim SMCs result in more SMCs, thereby an increase in total capacitance of the SMC layer. All else equal, a reduction in SMC width therefore slightly impairs electrical spread along the vessel wall and vice versa.

In Fig. 7, simulations show electrical spread across a single branch point initiated using current injection ð  900 pAÞ into ECs of one (black curves) or two vessels (gray curves). In Fig. 7A, the three vessels are identical except for the connectivity at the branch point (i.e. all other parameters are identical across simulations). In Fig. 7B, one vessel is shorter than the others but network lengths in Fig. 7A and B are equal. Our basic assumption is that ECs (the main electrical conduits) from each vessel connect to ECs from the other vessels in proportion to the number of ECs in each vessel (weighted division, see Fig. 4F). This type of connectivity, however, causes a substantially larger Vm -jump across the branch point compared to all-to-all coupling (compare full vs. dashed curve in Fig. 7A) despite the total number of EC:EC couplings at the branch point in this case is identical (see number above bars in Fig. 7C). Stimulation of two vessels (both  900 pA, gray curves) causes a stronger overall hyperpolarization of the entire branched network (compare gray to black curves) and the stronger stimulation only increases the Vm -jump. Shortening of the stimulated vessel (Fig. 7B) also causes a larger perturbation of Vm at the branch point and therefore produces a larger Vm -jump, illustrating the additional dependency on vessel morphology and stimulation strength. As summarized in Fig. 7C, all-to-all coupling at the branch point produces a smaller Vm -jump compared to coupling of vessels by division of ECs irrespective of the total number of couplings. Division of ECs at the branch point increases resistance between vessels as only some cells are directly coupled compared to all-to-all connectivity. This difference in resistance is observed as increased hyperpolarization of stimulated vessels in the simulations of weighted division connectivity of ECs (compared to dashed curves, all full curves are hyperpolarized in stimulated vessels and depolarized in unstimulated vessels, see Fig. 7A and B). Across the single branch point, the size of the Vm -jump may seem small (  5–10% of the total electrical decay). However, the effect of numerous branch points in a vascular network is additive, i.e. branching effects may be pronounced in the microcirculation.

3.2. Across a single branch point

3.3. At network level

From an electrical signaling perspective, a branch point represents a structural obstacle as current from one vessel is divided into other branches. The increased mass of cells at the branch point results in a larger overall capacitance, i.e. a larger electrical decay is expected across a branch point compared to the decay within a vessel (in the following the former is called a Vm -jump).

To illustrate electrical spread in bifurcating networks, an artificial tree-like network and an artificial arcading network were constructed as shown in Fig. 8. All vessel segments are identical and the networks differ only in topology, i.e. where the vessels couple to one another. Identical amounts of current ð  900 pAÞ were injected into either 3 upstream vessels (Fig. 8A and C),

8

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

stream stimulation. The propensity of a stimulated “capillary” to signal the feed-artery, controlling flow into the entire network, therefore depends on the number of branches in between and (to a lesser degree) on the total distance. Collectively, network topology influences (1) electrical sensitivity to stimulation, (2) the ratio between current across the plasma membrane (current dissipation) and current through GJs (dilution of charge within the network), and (3) the propensity of capillaries to signal the feed artery.

4. Discussion

Fig. 6. Single vessel level: cell size affects electrical communication through connectivity. Changing EC length, EC width or SMC width (A) gives rise to substantial changes in electrical conduction. Electrical spread was initiated via hyperpolarizing current injection (900 pA, B) or voltage clamp (  55 pA, C) of ECs of the first segment of a 2000 μm vessel. Conduction is summarized as a bar showing the distal Vm -response (Vm;distal  Vm;rest ) in black and the local Vm -response (Vm;local  Vm;rest ) in light gray. The percent-wise deviation from control of the distal Vm -response is shown within each bar. Cell size determines the number of cells in the circumferential and longitudinal directions which, in turn, partly determines the overall coupling pattern and hence, conduction.

3 neighboring downstream vessels (Fig. 8B and E), or 3 nonneighboring downstream vessels (Fig. 8C and F). As predicted, most electrical decay tends to occur at branch points (less apparent color change within individual vessels) in each simulation, i.e. due to division of ECs at the branch point. However, averaging across all vessels in the network, the stimulation regimes caused similar overall hyperpolarization in each network (  42:5 mV 70:12 mV with networks in the nonneighboring stimulation regimes being more hyperpolarized). Yet, stimulation of neighboring vessels that reside within subnetworks of low connectivity display less charge distribution across the network, i.e. such areas are more sensitive to electrical stimulation (compare Fig. 8B and E). Contrary, stimulation of nonneighboring vessels is readily diluted across the network (see Fig. 8C and F). Fig. 8C and F illustrates that fewer branches and short distances between downstream and upstream vessels produce more hyperpolarization of the feed-artery upon down-

A novel strategy for connectivity generation in computational models was presented and used to show that network topology as well as morphology and connectivity of ECs (locally and at branch points) delimits the electrical signaling capacity within vascular networks. Cell–cell layout and interactions may give rise to anisotropy at the local level and indirectly determine sensitivity of electrical perturbation at the global level. The implications of structure in extensive dynamical systems are usually hard to address both experimentally (as structure cannot be changed independently) and numerically as changes in connectivity entail model reorganization. The present generative strategy (Hald et al., 2013) imposes no restrictions on dynamics nor generative algorithms, i.e. the strategy can be applied to most biological systems. Here, a ‘vascular model generator’ was developed and used to construct (identical) dynamical models of electrical signaling in vascular networks of different size and connectivity. At a single vessel level, electrical conduction was enhanced by (a) increasing the number and longitudinal orientation of GJ couplings and (b) incorporating long and thin ECs. At a branch point, a weighted EC:EC coupling distribution between branches displays increased electrical resistance compared to the “all-to-all” coupling distribution which favors conduction into small branches. Finally, the network topology influences regional electrical sensitivity and thereby the effective propensity of a downstream signal to reach upstream vessel segments that effectively control overall network flow. As structural parameters delimits electrical communication, structure may therefore be expected to feedback on local electrical behavior, i.e. the local topology may partially shape the dynamical properties of a vascular bed. 4.1. The generative modeling approach as applied to the vasculature The generative approach to connectivity embodied in the vascular model generator extends our previously published strategy (Hald et al., 2013) that promotes separation of dynamics within and between biological entities, i.e. internal dynamics and dynamics of interactions between entities are described individually. This allows for arbitrary assignment of connections and construction of larger hierarchical structures (see Fig. 2). In terms of the vasculature, connectivity generation by algorithm(s) can be summarized in the following: (1) Within a vessel segment: if the morphologies and layout of the vessel and constituent vascular cells are known, the connectivity within a vessel can be calculated. (2) Between a vessel segments: if the connectivity of vessel segments within a vascular network is known (graph generated either via reconstruction or synthetic), the particular connectivity of ECs and SMCs at a branch point can be calculated. Consequently, arbitrary vascular networks can be generated depending on input. This underlines one benefit of developing a model generator program instead of building a dedicated program for a specific single model. Here, the vascular model generator used a directed graph to handle vessel connectivity (as described by the user) whereas cell–cell connectivity was generated. To produce a

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

9

Fig. 7. Electrical spread is sensitive to connectivity at branch points. In all simulations shown, 900 pA current was injected into ECs of the first segment of a single vessel (black curves) or two vessels (gray curves), see insets, giving rise to electrical spread across a branch point to neighboring vessels. At the branch point, ECs from each vessel were divided by “weighted division” (full curve), i.e. ECs of each partition have 6 couplings to ECs from a single partition of each other vessels, or were coupled “all-to-all” (dashed curve), i.e. each EC couples to one EC from every other vessels (see Fig. 4F). (A) Three identical vessels were connected in a single branch point. (B) One vessel is shorter than the others. However, the cumulative network length is the same as in A. (C) Summary. The Vm -decay across the branch point (Vm -jump) as a function of connectivity resulting from stimulation of a single vessel (dark gray) or two vessels (light gray). Below each bar, “all” indicates all-to-all connectivity whereas a number indicates the number of EC:EC couplings using weighted division connectivity. The number above each bar denotes the total number of couplings in the branch point. Stronger stimulation (dual vs single or unequal vs equal vessel size) causes larger Vm -jump. Division of ECs at the branch point also causes larger Vm -jumps compared to allto-all coupling – even if more couplings are present at the branch point.

b

c

Tree network

a

Upstream stimulation

e

Downstream stimulation

f

Arcading network

d

Downstream stimulation

200 µm

Fig. 8. Electrical spread in two artificial networks. Using 61 vessels with identical properties, a tree-like network (upper panel) and a planar arcading network (lower panel) were constructed.  900 pA current injections into ECs of the first circumferential segment of 3 selected vessels (see red arrow-heads) produced different patterns of electrical spread within the network. Red color indicates vessels at Vm;rest , i.e. unperturbed by the stimulation. Hyperpolarization is indicated by increasingly blue color. (A) and (D) Current injection in 3 upstream vessels. (B) and (E) Current injection in 3 downstream neighboring vessels. (C) and (F) Current injection in 3 upstream nonneighboring vessels. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

network model, the minimal input to the vascular model generator is reduced to the node coordinates (in 2D or 3D) delineating vessels. Resting diameter (and length), number of cell layers, and cell dimensions are optionally assigned to each vessel segment (see Python routines in the Supplemental Information). The connectivity generating algorithm(s) can be made more or less complicated/accurate depending on experimental knowledge or the user-defined problem without invalidating the central ideas. 4.2. Structure delimit electrical communication in vascular networks The utility of the vascular model generator was illustrated by generating series of structurally different systems while keeping electrical kinetics within and between vascular cells constant. Overall gap junction conductance is an extensive property: In the vascular conduction field, model developers make more or less arbitrary choices regarding GJ couplings and morphology. Models of conduction either use a constant value for individual GJ conductance (Diep et al., 2005; Wolfle et al., 2011; Kapela et al., 2010) or employ GJ conductances that scale with the membrane surface

area (Crane et al., 2001; Jacobsen et al., 2007; Koenigsberger et al., 2004). As illustrated in Figs. 5 and 6, electrical conduction along a single vessel is not simply a function of the mean conductance of a gap junction plaque. Rather, electrical conduction also depends on the average number of couplings per cell (Fig. 5B), the predominant orientation of GJs (anisotropy in connectivity, Fig. 5C), and the physical dimensions of vascular cells (anisotropy based on structure, Fig. 6), particularly ECs. In conduction models, these parameters are vaguely defined, obscuring the relation between the chosen GJ conductance and electrical spread. Whereas the physical dimensions of vascular cells can be relatively easily measured (see below), experimental evidence regarding the predominant orientation and average number of couplings per vascular cell is lacking. Scaling with number of couplings dissociates simulations of conduction from this value (Fig. 5B). However, introduction of anisotropy in connectivity can still influence conduction (Fig. 5C): as more couplings are introduced along the longitudinal axis, the relative coupling strength increases longitudinally but decreases circumferentially. Equivalently, if e.g. longitudinally oriented cells display extensive overlap, coupling

10

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

strength (although uniformly distributed within an individual cell) also increases in the longitudinal direction. To the best of our knowledge, GJs are not known to show preferential locations within individual vascular cells. However, both SMCs and ECs overlap extensively (Haas and Duling, 1997; Behringer et al., 2012; Martinez-Lemus et al., 2004), i.e. some degree of structural anisotropy is probable. To emulate overlapping cells, we chose 6 coupling neighbors per cell as a standard (note that coupling pattern is not equivalent to the physical position of cells). Scaling requires often-used experimental estimates of GJ conductances (Yamamoto et al., 2001; Lidington et al., 2000) to be divided by the number of couplings. Therefore, the scaled GJ conductance produced a strong conduction decay incommensurate with an exponential decay constant around 1–2 mm that is routinely observed experimentally (compare dashed curve to gray shaded area in Fig. 5C) (Xia and Duling, 1995; Segal et al., 1989; Emerson et al., 2002). Decay constant above  1 mm was obtained by multiplying GJ conductances per cell by 8 (solid curves in Fig. 5C), suggesting that experimental estimates of GJ conductance are somewhat low. As experimental measures of IV-curves of the SMC and EC are relatively well-defined (see Fig. 3A), we propose that average coupling resistances between cells (Yamamoto et al., 2001; Lidington et al., 2000) are better estimated using a proper electrical model and fitting average GJ conductance per cell to measure decay along the vascular wall. Electrical spread in vessel segments is affected by cell morphology: In contrast to the SMC layer, the endothelium is known to be an efficient electrical pathway due to strong coupling conductance (Diep et al., 2005; Hald et al., 2012) and a parallel cell orientation to the vessel length axis (Haas and Duling, 1997; Behringer et al., 2012; Hill, 2012; Diep et al., 2005; Hald et al., 2012), i.e. as electrical charge spreads along the vascular wall it faces the lowest possible number of GJs (cell-cell couplings) with the highest conductance in the EC layer. The same argument explains why increasing average EC length enhances conduction (Fig. 6). Reducing EC width (slim ECs) increases the number of circumferential ECs as well as the total membrane area between neighboring ECs. As this area is assumed to be proportional to GJ conductance per cell, a reduction in EC width corresponds to an increase in overall coupling strength and hence, conduction efficacy (Fig. 6B). Finally, slim cells increase the probability of more cells being stimulated by some locally acting external factor and local Vm -response is proportional to the number of stimulated cells (Fig. 6A). As the SMC layer is ill-suited for longitudinal electrical spread, the effect of changing average SMC size is limited. However, a decrease in SMC width increases in the number of SMCs in the vessel segment. As total MEGJ conductance is independent of the number of SMCs (same surface area between layers), more SMCs correspond to a higher total capacitance of the SMC layer. A higher capacitive load results in larger conduction decay (Fig. 6) (Kapela et al., 2010). Despite the idealized nature of the computational model: assignment of identical cell dimensions and kinetic properties to every EC and SMC, the effect of varying “average” cell size on conduction decay is evident. In the literature, EC length and width is observed to vary somewhat within a vascular network but considerably between different networks (Haas and Duling, 1997; Sandow et al., 2003a). Also, cell size changes during development (Sandow et al., 2004). Finally, SMC width has been reported to vary between 2 μm and 10 μm in different networks (Haas and Duling, 1997; Gattone et al., 1986; Sandow et al., 2003b). With all else equal, this implies that vascular beds by their structure alone are endowed with different structural dispositions for effective electrical communication. In conclusion, electrical spread along a single vessel is significantly influenced by the particular connectivity and/or morphology and is not only determined by GJ and

plasma membrane conductances of individual cells (Hald et al., 2015; Diep et al., 2005). Connectivity influences electrical spread across a branch point: Compared to single vessels, electrical spread is known to increase decay across branch points (Kurjiaka and Segal, 1995; Segal and Neild, 1996). This is reproduced by the simulations (see Fig. 7). A branch point is equivalent to a discontinuity in the total mass and capacitance (the mass of a parent vessel is not conserved in the mass of the two or more daughter vessels). Thus, the electrical conduction profile across a branch point must invariably show a larger Vm -decay (Vm -jump) as compared to the incremental decay observed between cells within a vessel segment. With ECs being the primary electrical conduits, charge spread across a branch point depends critically on the local EC:EC connectivity (Fig. 7). Naturally, a lower stimulation strength reduces the Vm -jump. The Vm -jump also decreases with more couplings but notably also with better cross-coupling between vessels. Despite very limited data from the literature, some level of cross-coupling is observed in the SMC layer (Ushiwata and Ushiki, 1990). As ECs are oriented along the vessel, branching must incur some distortion of ECs, suggesting some cross-coupling ECs. On the other hand, branching between large and small arterioles limits the extent of crosscoupling, preventing too much current flow into the small vessel. That is, extensive cross-coupling of ECs between branches corresponds to a disproportional charge spread into small vessels. The subtle effect of connectivity across a branch point becomes important in vascular networks as the branch point precisely affects charge distribution among vessels. Indeed, in Fig. 8, the effect of branch points is observed as points where most electrical decay appears. Electrical spread in vascular networks: Due to the lack of freely available databases, artificial networks were constructed to illustrate the dependency of electrical spread on network topology. The two networks in Fig. 8 are composed of identical vessel segments and show that networks with a high degree of “connectedness” readily distribute current within the network at the expense of producing weak local perturbations. Conversely, vascular networks with low connectivity become sensitive to electrical stimuli. Thus, the Vm of some regions in a network may be highly sensitive to electrical stimuli. Arcading networks and/or networks with many collaterals are typically highly “connected”, while bifurcating tree-like networks tend to have low connectivity (Hirsch, 2012). As IV-curves are often monotonically increasing functions (Fig. 3A), a stronger Vm -perturbation away from Vm;rest is therefore typically associated with higher current dissipation (particularly upon depolarizing stimuli). Consequently, upon identical initial stimulation, a vascular region with low connectivity tends to produce a strong local Vm -response that shows large current dissipation, whereas stimulation of a highly interconnected region tends to preserve charge better within the system at the expense of diminishing the local electrical signal. In excess, both scenarios are detrimental to effective electrical communication in vascular networks. In contrast to the artificial networks considered in Fig. 8, differences in electrical spread may be more pronounced in real networks where structure and morphology changes from capillaries to feed arteries (i.e. across Strahler orders). Any network reconstruction methodology (outside the scope of this paper, e.g. Cassot et al., 2006; Reichold et al., 2009; Blinder et al., 2013; Kim et al., 2012), however, could be used as input to the vascular model generator to develop dynamical models based on real morphometric network data (Fig. S1 in the Supplemental Information illustrates a cerebral network generated from a small freely available data set). Automatic construction of dynamical models from different vascular beds could be used to study how structure limits the flow of information within the microvasculature. For

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

example, electrical spread arising from stimulation of different regions may be disproportionally dispositioned to reach the feed artery (compare e.g. Fig. 8B and C) that gates blood flow into the network. Such comparisons, however, require freely available data-sets which are currently scarce despite the abundance of reconstruction techniques of 3D networks in the literature. 4.3. Model limitations and future directions Simple dynamical models of electrical behavior in ECs and SMCs were used due to the focus on connectivity generation. These models describe the aggregate current of every ion channel/ exchanger/pump in the plasma membrane as a single 1dimensional IV-curve. This is obviously a gross simplification of the cumulative behavior of ion transporters and neglect most important vasomotor responses, signaling cascades, and feedback mechanisms but it has the following virtues: (1) IV-curves are very simple to describe and use. They derive directly from experimental patch-clamp experiments using voltage-ramps in physiological solutions (other solutions, however, might change the Nernst equilibrium and consequently the IV-curve). Whereas the influence of any chemical or mechanical stimulation is not directly included, such stimuli are to some extent contained in the IVcurve insofar that these stimuli are “autocrine” and electrical. An IV-curve describes current across the membrane as a function of voltage (here, the only variable in the model). However, the IVcurves of ECs and SMCs located in different regions of the network may be different in terms of electrical behavior, but this is a subject for future studies. Here, the effect of structure on electrical spread in structurally different vascular networks was studied and simple IV-models suffice. An ohmic representation of GJ currents is normally used. It is well-known that both hyperpolarization and depolarization can spread between coupled (vascular) cells but if the voltagedifference becomes too large, the GJ will close (the exact behavior depends on connexin sub-type expression). Within a “physiologically” relevant region, i.e. ΔVm 7 15 mV, the normalized GJ conductance should be approximately at unity with minimal anisotropy (Vogel and Weingart, 2002). A set of structural limitations also applies. In general, the strategy describes spatial location indirectly through the connectivity of the individual cell. The accuracy of the generator algorithms is of course crucial for a representative model, but their influence can be directly assessed by comparison either to experimental structures or between different generator algorithms. In most cases, identical connectivity patterns will be applied to individual cells but this can be ameliorated somewhat through the use of random generators, see Hald et al. (2013). Likewise, the parameters designating volume, surface area or kinetic parameters can be subject to a defined degree of heterogeneity (not used here). As the vascular model generator accepts simple morphometric data, more detailed network descriptions of real vascular networks from different tissues should enable a better understanding of electrical flow patterns. Finally, the dynamics included can easily be extended e.g. by inclusion of wall diameter and/or blood flow models. Thus, the model generator can be a powerful research tool to study flow patterns in vascular networks.

5. Conclusions The present paper described a novel strategy to generate connectivity in models of structurally extensive biological systems, presented versatile software that employ the strategy in dynamical models of electrical signaling in vascular networks, and

11

showed that electrical spread in vascular networks is shaped by connectivity and morphology. Use of arbitrary pattern-generators of connectivity and/or ODE sets is possible via separation of dynamics and connectivity. Here, simulations of series of structurally different models show that electrical conduction is enhanced by increasing the number of longitudinally oriented couplings and/or in vessels harboring long and thin ECs (and/or thick SMCs). Second, branch point connectivity and network topology determine regional sensitivity to electrical stimuli and, in turn, the propensity of regional stimuli to elicit perturbation of the feed artery, gating the overall network flow. The software could be useful to study altered information flows upon the structural changes associated with vascular disease, e.g. vascular remodeling and/or rarefaction observed in hypertension or as a template for more sophisticated models of either vascular dynamics or similar biological systems, e.g. the biliary system or lymphatics.

Acknowledgments BOH is supported by the Danish Council for Independent Research (Sapere Aude program: DFF – 1333-00172 and DFF – 1331-00731B) and the Dynamical Systems Interdisciplinary Network, University of Copenhagen. The author would like to thank P. G. Sørensen, J.C. Brasen, and J.C.B. Jacobsen for helpful discussions. Author Contributions: BOH designed and performed the research and wrote the paper.

Appendix A. Supplementary data Supplementary data associated with this paper can be found in the online version at http://dx.doi.org/10.1016/j.jtbi.2016.03.032.

References Behringer, E.J., Socha, M.J., Polo-Parada, L., Segal, S.S., 2012. Electrical conduction along endothelial cell tubes from mouse feed arteries: confounding actions of glycyrrhetinic acid derivatives. Br. J. Pharmacol. 166 (2), 774–787. http://dx.doi. org/10.1111/j.1476-5381.2011.01814.x. Blinder, P., Tsai, P.S., Kaufhold, J.P., Knutsen, P.M., Suhl, H., Kleinfeld, D., 2013. The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow. Nat. Neurosci. 16 (7), 889–897. http://dx.doi.org/10.1038/ nn.3426. Cassot, F., Lauwers, F., Fouard, C., Prohaska, S., Lauwers-Cances, V., 2006. A novel three-dimensional computer-assisted method for a quantitative study of microvascular networks of the human cerebral cortex. Microcirculation 13 (1), 1–18. http://dx.doi.org/10.1080/10739680500383407. Cole, W.C., Welsh, D.G., 2011. Role of myosin light chain kinase and myosin light chain phosphatase in the resistance arterial myogenic response to intravascular pressure. Arch. Biochem. Biophys. 510 (2), 160–173. http://dx.doi.org/10.1016/j. abb.2011.02.024. Cornelis, H., Rodriguez, A.L., Coop, A.D., Bower, J.M., 2012. Python as a federation tool for genesis 3.0. PLoS One 7 (1), e29018. http://dx.doi.org/10.1371/journal. pone.0029018. Crane, G.J., Hines, M.L., Neild, T.O., 2001. Simulating the spread of membrane potential changes in arteriolar networks. Microcirculation 8 (1), 33–43. Diep, H.K., Vigmond, E.J., Segal, S.S., Welsh, D.G., 2005. Defining electrical communication in skeletal muscle resistance arteries: a computational approach. J. Physiol. 568 (Pt 1), 267–281. http://dx.doi.org/10.1113/jphysiol.2005.090233. Emerson, G.G., Neild, T.O., Segal, S.S., 2002. Conduction of hyperpolarization along hamster feed arteries: augmentation by acetylcholine. Am. J. Physiol. Heart Circ. Physiol. 283 (1), H102–H109. http://dx.doi.org/10.1152/ajpheart.00038.2002. Gattone, V.H., Miller, B.G., Evan, A.P., 1986. Microvascular smooth muscle cell quantitation from scanning electron microscopic preparations. Anat. Rec. 216 (3), 443–447. http://dx.doi.org/10.1002/ar.1092160315. Haas, T.L., Duling, B.R., 1997. Morphology favors an endothelial cell pathway for longitudinal conduction within arterioles. Microvasc. Res. 53 (2), 113–120. http: //dx.doi.org/10.1006/mvre.1996.1999. Hald, B.O., Welsh, D.G., Holstein-Rathlou, N.-H., Jacobsen, J.C.B., 2015. Origins of variation in conducted vasomotor responses. Pflugers Arch 467 (10), 2055–2067.

12

B.O. Hald / Journal of Theoretical Biology 399 (2016) 1–12

Hald, B.O., Smrcinova, M., Sørensen, P.G., 2012. Influence of cyanide on diauxic oscillations in yeast. FEBS J. 279 (23), 4410–4420. http://dx.doi.org/10.1111/ febs.12030. Hald, B.O., Jensen, L.J., Sørensen, P.G., Holstein-Rathlou, N.-H., Jacobsen, J.C.B., 2012. Applicability of cable theory to vascular conducted responses. Biophys. J. 102 (6), 1352–1362. http://dx.doi.org/10.1016/j.bpj.2012.01.055. Hald, B.O., Hendriksen, M.G., Sørensen, P.G., 2013. Programming strategy for efficient modeling of dynamics in a population of heterogeneous cells. Bioinformatics 29 (10), 1292–1298. http://dx.doi.org/10.1093/bioinformatics/btt132. Hald, B.O., Jacobsen, J.C.B., Sandow, S.L., Holstein-Rathlou, N.-H., Welsh, D.G., 2014. Less is more: minimal expression of myoendothelial gap junctions optimizes cell-cell communication in virtual arterioles. J. Physiol. 592 (Pt 15), 3243–3255. http://dx.doi.org/10.1113/jphysiol.2014.272815. Hill, M.A., Davis, M.J., Meininger, G.A., Potocnik, S.J., Murphy, T.V., 2006. Arteriolar myogenic signalling mechanisms: implications for local vascular function. Clin. Hemorheol. Microcirc. 34 (1–2), 67–79. Hill, C.E., 2012. Long distance conduction of vasodilation: a passive or regenerative process?. Microcirculation 19 (5), 379–390. http://dx.doi.org/10.1111/ j.1549-8719.2012.00169.x. Hirsch, S., Reichold, J., Schneider, M., Székely, G., Weber, B., 2012. Topology and hemodynamics of the cortical cerebrovascular system. J. Cereb. Blood Flow Metab. 32 (6), 952–967. http://dx.doi.org/10.1038/jcbfm.2012.39. Jackson, W.F., Huebner, J.M., Rusch, N.J., 1997. Enzymatic isolation and characterization of single vascular smooth muscle cells from cremasteric arterioles. Microcirculation 4 (1), 35–50. Jacobsen, J.C.B., Gustafsson, F., Holstein-Rathlou, N.-H., 2003. A model of physical factors in the structural adaptation of microvascular networks in normotension and hypertension. Physiol. Meas. 24 (4), 891–912. Jacobsen, J.C.B., Aalkjaer, C., Nilsson, H., Matchkov, V.V., Freiberg, J., HolsteinRathlou, N.-H., 2007. A model of smooth muscle cell synchronization in the arterial wall. Am. J .Physiol. Heart Circ. Physiol. 293 (1), H229–H237. http://dx. doi.org/10.1152/ajpheart.00727.2006. Kapela, A., Nagaraja, S., Tsoukias, N.M., 2010. A mathematical model of vasoreactivity in rat mesenteric arterioles. ii. conducted vasoreactivity. Am. J. Physiol. Heart Circ. Physiol. 298 (1), H52–H65. http://dx.doi.org/10.1152/ ajpheart.00546.2009. Kim, E., Stamatelos, S., Cebulla, J., Bhujwalla, Z.M., Popel, A.S., Pathak, A.P., 2012. Multiscale imaging and computational modeling of blood flow in the tumor vasculature. Ann. Biomed. Eng. 40 (11), 2425–2441. http://dx.doi.org/10.1007/ s10439-012-0585-5. Koenigsberger, M., Sauser, R., Lamboley, M., Bény, J.-L., Meister, J.-J., 2004. Ca2 þ dynamics in a population of smooth muscle cells: modeling the recruitment and synchronization. Biophys. J. 87 (1), 92–104. http://dx.doi.org/10.1529/ biophysj.103.037853. Kurjiaka, D.T., Segal, S.S., 1995. Conducted vasodilation elevates flow in arteriole networks of hamster striated muscle. Am. J. Physiol. 269 (5 (Pt 2)), H1723–H1728. Lidington, D., Ouellette, Y., Tyml, K., 2000. Endotoxin increases intercellular resistance in microvascular endothelial cells by a tyrosine kinase pathway. J. Cell Physiol. 185 (1), 117–125. http://dx.doi.org/10.1002/1097-4652(200010) 185:13.0.CO;2-7. Longden, T.A., Hill-Eubanks, D.C., Nelson, M.T., 2016. Ion channel networks in the control of cerebral blood flow. J. Cereb. Blood Flow Metab. 36 (3), 492–512. Marsh, D.J., Wexler, A.S., Brazhe, A., Postnov, D.E., Sosnovtseva, O.V., HolsteinRathlou, N.-H., 2013. Multinephron dynamics on the renal vascular network. Am. J. Physiol. Ren. Physiol. 304 (1), F88–F102. http://dx.doi.org/10.1152/ ajprenal.00237.2012. Martinez-Lemus, L.A., Hill, M.A., Bolz, S.S., Pohl, U., Meininger, G.A., 2004. Acute mechanoadaptation of vascular smooth muscle cells in response to continuous arteriolar vasoconstriction: implications for functional remodeling. FASEB J. 18 (6), 708–710. http://dx.doi.org/10.1096/fj.03-0634fje. Mehta, D., Malik, A.B., 2006. Signaling mechanisms regulating endothelial permeability. Physiol. Rev. 86 (1), 279–367. http://dx.doi.org/10.1152/ physrev.00012.2005.

Nilius, B., Droogmans, G., 2001. Ion channels and their functional role in vascular endothelium. Physiol. Rev. 81 (4), 1415–1459. Pries, A.R., Secomb, T.W., 2014. Making microvascular networks work: angiogenesis, remodeling, and pruning. Physiology 29 (6), 446–455. http://dx.doi.org/ 10.1152/physiol.00012.2014. Reichold, J., Stampanoni, M., Keller, A.L., Buck, A., Jenny, P., Weber, B., 2009. Vascular graph model to simulate the cerebral blood flow in realistic vascular networks. J. Cereb. Blood Flow Metab. 29 (8), 1429–1443. http://dx.doi.org/10.1038/ jcbfm.2009.58. Rieger, H., Welter, M., 2015. Integrative models of vascular remodeling during tumor growth. Interdiscip. Rev. Syst. Biol. Med. 7 (3), 113–129. http://dx.doi.org/ 10.1002/wsbm.1295. Rijen, H.V., van Kempen, M.J., Analbers, L.J., Rook, M.B., van Ginneken, A.C., Gros, D., Jongsma, H.J., 1997. Gap junctions in human umbilical cord endothelial cells contain multiple connexins. Am. J. Physiol. 272 (1 (Pt 1)), C117–C130. Sandow, S.L., Looft-Wilson, R., Doran, B., Grayson, T.H., Segal, S.S., Hill, C.E., 2003a. Expression of homocellular and heterocellular gap junctions in hamster arterioles and feed arteries. Cardiovasc. Res. 60 (3), 643–653. Sandow, S.L., Bramich, N.J., Bandi, H.P., Rummery, N.M., Hill, C.E., 2003b. Structure, function, and endothelium-derived hyperpolarizing factor in the caudal artery of the shr and wky rat. Arterioscler. Thromb. Vasc. Biol. 23 (5), 822–828. http: //dx.doi.org/10.1161/01.ATV.0000067425.06019.D7. Sandow, S.L., Goto, K., Rummery, N.M., Hill, C.E., 2004. Developmental changes in myoendothelial gap junction mediated vasodilator activity in the rat saphenous artery. J. Physiol. 556 (Pt 3), 875–886. http://dx.doi.org/10.1113/ jphysiol.2003.058669. Schneider, M., Hirsch, S., Weber, B., Székely, G., Menze, B.H., 2014. Tgif: topological gap in-fill for vascular networks-a generative physiological modeling approach. Med. Image Comput. Comput. Assist. Interv. 17 (Pt 2), 89–96. Schnitzler, M.M.y., Derst, C., Daut, J., Preisig-Müller, R., 2000. Atp-sensitive potassium channels in capillaries isolated from guinea-pig heart. J. Physiol. 525 (Pt 2), 307–317. Segal, S.S., Neild, T.O., 1996. Conducted depolarization in arteriole networks of the guinea-pig small intestine: effect of branching of signal dissipation. J. Physiol. 496 (Pt 1), 229–244. Segal, S.S., Damon, D.N., Duling, B.R., 1989. Propagation of vasomotor responses coordinates arteriolar resistances. Am. J. Physiol. 256 (3 Pt 2), H832–H837. Segal, S.S., 2005. Regulation of blood flow in the microcirculation. Microcirculation 12 (1), 33–45. http://dx.doi.org/10.1080/10739680590895028. Somlyo, A.P., Wu, X., Walker, L.A., Somlyo, A.V., 1999. Pharmacomechanical coupling: the role of calcium, g-proteins, kinases and phosphatases. Rev. Physiol. Biochem. Pharmacol. 134, 201–234. Stoma, S., Fröhlich, M., Gerber, S., Klipp, E., 2011. Stse: spatio-temporal simulation environment dedicated to biology. BMC Bioinform. 12, 126. http://dx.doi.org/ 10.1186/1471-2105-12-126. Tran, C.H.T., Vigmond, E.J., Goldman, D., Plane, F., Welsh, D.G., 2012. Electrical communication in branching arterial networks. Am. J. Physiol. Heart Circ. Physiol. 303 (6), H680–H692. http://dx.doi.org/10.1152/ajpheart.00261.2012. Ushiwata, I., Ushiki, T., 1990. Cytoarchitecture of the smooth muscles and pericytes of rat cerebral blood vessels. a scanning electron microscopic study. J. Neurosurg. 73 (1), 82–90. http://dx.doi.org/10.3171/jns.1990.73.1.0082. Vogel, R., Weingart, R., 2002. The electrophysiology of gap junctions and gap junction channels and their mathematical modelling. Biol. Cell 94 (78), 501–510. Wölfle, S.E., Chaston, D.J., Goto, K., Sandow, S.L., Edwards, F.R., Hill, C.E., 2011. Nonlinear relationship between hyperpolarisation and relaxation enables long distance propagation of vasodilatation. J. Physiol. 589 (Pt 10), 2607–2623. http: //dx.doi.org/10.1113/jphysiol.2010.202580. Xia, J., Duling, B.R., 1995. Electromechanical coupling and the conducted vasomotor response. Am. J. Physiol. 269 (6 (Pt 2)), H2022–H2030. Yamamoto, Y., Klemm, M.F., Edwards, F.R., Suzuki, H., 2001. Intercellular electrical communication among smooth muscle and endothelial cells in guinea-pig mesenteric arterioles. J. Physiol. 535 (Pt 1), 181–195.