Network architecture associated with the highly specialized hindlimb of frogs

Network analyses have been increasingly used in the context of comparative vertebrate morphology. The structural units of the vertebrate body are treated as discrete elements (nodes) of a network, whose interactions at their physical contacts (links) determine the phenotypic modules. Here, we use the network approach to study the organization of the locomotor system underlying the hindlimb of frogs. Nodes correspond to fibrous knots, skeletal and muscular units. Edges encode the ligamentous and monoaxial tendinous connections in addition to joints. Our main hypotheses are that: (1) the higher centrality scores (measured as betweenness) are recorded for fibrous elements belonging to the connective system, (2) the organization of the musculoskeletal network belongs to a non-trivial modular architecture and (3) the modules in the hindlimb reflect functional and/or developmental constraints. We confirm all our hypotheses except for the first one, since bones overpass the fibrous knots in terms of centrality. Functionally, there is a correlation between the proximal-to-distal succession of modules and the progressive recruitment of elements involved with the motion of joints during jumping. From a developmental perspective, there is a correspondence between the order of the betweenness scores and the ontogenetic chronology of hindlimbs in tetrapods. Modular architecture seems to be a successful organization, providing of the building blocks on which evolution forges the many different functional specializations that organisms exploit.


Introduction
Network analyses have been increasingly used in the context of comparative vertebrate morphology due to their undoubtedly fresh perspective regarding issues such as anatomical topology, morphogenetic integration, and modularity [1,2]. Addressing anatomical systems as networks and exploring their patterns of connections can lead to useful insights. For instance, one of the most interesting properties of networks is intermediacy in information flow, and the betweenness score as a measure of centrality can be used for detecting the bridging role a given node can achieve in the flow of information through the network [3]. The development of a novel methodology for anatomical studies, such as anatomical network analysis (AnNA, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 system, (2) the musculoskeletal network exhibits a modular organization or non-trivial architecture, (3) the modules in the hindlimb reflect functional and/or developmental constraints, and these restrictions are not necessarily overlapped.
Networks are collections of discrete entities or elements (vertices or nodes) connected by some relationship of interest (links or edges). In our case, nodes correspond to anatomical elements (bones, muscles and fibrous knots) and edges encode physical contacts among such elements (e.g. tendinous and ligamentous junctions). Fig 1 shows the actual anatomical substrate from where the network is constructed.
Fibrous knots are here defined as critical portions of connective tissues from which multiple (more than two) sets of identically aligned collagen fibers seem to branch out (Fig 2). They are hypothesized to act as gateways in the transmission of mechanical signals. Consequently, edges incident on the fibrous knots represent tendinous segments, or equivalently, they account for fibrous substrates exposed to uniaxial mechanical stress in the actual anatomical system. The adjacency list, with the respective set of neighbors of each node, can be found in Table 2.
The anatomical elements of this study include intrinsic components of the pelvic girdle, stylopodium, and zeugopodyum, and extrinsic elements of the autopodium; except for the short muscles of the acetabular joint (quadratus femoris, gemelus, and obturador internus) (Fig 1). Due to the fact that the distal iliac area, ischium and pubis emerge as a single mesenchymal matrix during development [62,63], and are fused into one element in the adult stage, they were merged into the 'acetabulum' element of the matrix. We consider the same point of origin for all the flexor tendons of the foot, because they integrate the same superficial tendinous layer (i.e. the aponeurosis plantaris) despite their double origin (flexor tendons of digits I, II, and part of III arise from the aponeurosis plantaris and those of digits III, IV and V from m. flexor brevis superficialis). This study is performed with specimens that were already housed in herpetological collections; thus, permissions to collect specimens or the following of ethical regulations were not necessary.

Network analysis
General concepts. A network is a collection of nodes (or vertices) and links (or edges) modeling the elements of a system and their relationships, respectively. Mathematically, networks are studied by the graph theory. Below, we introduce certain technical concepts in order to facilitate the lecture of the manuscript by a broad readership.
Links can be directed or undirected depending on the oriented nature of the relationship, and weighted or unweighted if the strength of the relationship is valued or not, respectively. A path in a graph represents an alternate sequence of vertices and edges from a certain origin to a certain destination by traversing edges. The first vertex of the first edge of a path is the origin and the second vertex of the last edge is the destination. Both origin and destination are called the endpoints of the path. The length of a path is the number of edges that it uses. Given a graph G, the distance d (x, y) between two vertices x and y is the length of the shortest (or geodesic) path from x to y, considering all the possible paths from x to y in G.
Betweenness. Centrality measures capture the relevance of the position of the individual nodes in the network. One of these quantitative indices is the betweenness score. Betweenness accounts for the frequency of occurrence of a given node in shortest paths between any pair of nodes in the network. It is useful for detecting the bridging role a given node can exhibit in the flow of information across the structure of the network. As nodes with high betweenness are removed, the remaining elements of the network are prone to be disconnected. In the context of muscle-skeletal networks, centrality indicates high exposure of elements to the propagation of force.
Modularity and network layout. Modularity refers to densely connected clusters of nodes embedded into a more global network. It quantifies the extent, relative to a null model network, to which vertices are grouped into dense sub-graphs [3]. The critical issue for quantifying it relies on the definition of a modularity matrix B = A-P, in which A is the adjacency matrix of the network that defines a graph with n vertices and m edges, and the matrix P contains the probability of two nodes to be connected by an edge under the prescriptions of a random configuration model. In other words, each element of P = [p ij ] corresponds to the probability of existence of an edge between vertices i and j in a random network in which the degrees of all vertices are the same as in the input graph. For unimodal networks, the usual null model sets such probability as proportional to the product of degrees k (number of actual neighbors) of the nodes involved, namely p ij = k i k j /(2m). A measure of modularity, called coefficient Q, is devised in [3] and consists of adding the entries of the modularity matrix for pairs of vertices belonging to the same group of a given vertex set partition. Thus, the modularity Q, for a given assignment g of vertices to groups or modules, reflects the extent to which edges are formed within modules instead of between modules, relative to a null model. Formally, modularity is defined as where the right hand factor is the Kronecker delta, which is equal to 1 if nodes i and j are in the same module and is otherwise 0. The partition, or division of vertices into modules, that maximizes Q is preferred. Since modularity optimization is NP-hard, heuristic procedures have been developed. Among the many available algorithms implemented in the igraph package of R, we achieve the best performance with the method of Brandes et al. [64] coded in a function named cluster_optimal. This method uses an Integer Linear Programming formulation to maximize modularity.   [60]. The Achilles tendon originates from the flexor digitorum communis (fdc) and then spreads out along the underside of the foot to form the aponeurosis plantaris, which in turn sends tendinous strands into each of the several digits. For this particular case, musculoskeletal connections between fdc and digits can be modeled through two alternatives, one more naïve (C) and one more realistic (D) network representation. In the naïve case, several edges are incident on the same muscular node suggesting that several tendons originate from such muscle. In the realistic case, just a single edge is incident on the muscular node, better reflecting the existence of the single Achilles tendon originated from the muscle under consideration.  Information about the topology of network relationships is encoded in the adjacency matrix, but structural patterns, such as modularity, can be easily revealed through appropriate network drawings. Force-directed algorithms for the placement of nodes within a rectangular frame of simple undirected graphs are noteworthy. Besides aesthetically pleasing, their output frequently reveals aspects about the organization of links, because such algorithms calculate the layout of a graph using only information contained within the structure of the graph itself. The layout method of Fruchterman & Reingold [65] belongs to that category of graph drawing and is inspired in spring forces proper of natural systems. In this method, the final location of a node is dictated by the interplay between repulsive and attractive forces between nodes. Repulsive forces act among all pairs of nodes, while attractive forces act between adjacent nodes. The final output aims to capture as much symmetry as possible and pursues to evenly distribute the vertices in the frame. Since vertices connected by an edge tend to be drawn near each other, modules can be recognized as proximate subsets of connected nodes. When dealing with sparse small graphs (such as ours), this method is highly recommended for cluster detection. Different analyses and graphics were performed through the add-on package igraph, available in R software [66].

Overall network layout
The network associated with the musculoskeletal system of the anuran hindlimb is showed in Fig 3. The two-dimensional arrangement of nodes aims to reflect the modular structure underlying the network. This layout, automatically obtained after applying the Fruchterman-Reinglod algorithm, largely mirrors the overall architecture of the hindlimb: the modular organization and proximal-distal arrangement of elements.
We compiled 27 muscles, 20 bones and 7 fibrous knots as nodes. With the exception of the cruralis and the extensor iliotibialis B, all pairs of muscular elements remain without a direct link between them. Additionally, 22 out of 27 muscles are of degree 2 reflecting origin and insertion into skeletal elements in a simple way.

Centrality
The system includes a few nodes with higher centrality values. Nine out of 54 nodes (about 17%) show the highest values of betweenness in the data set (larger than 100, Fig 4). The decreasingly ordered sequence of betweenness scores evokes an exponential decay. The osteological elements account for the majority of elements belonging to the set with higher centrality values: femur, tibiofibula, acetabulum, iliac shaft and tibiale. Comparatively, muscles yield the lower centrality values. All the elements related to the autopodium show low betweenness

Modularity
The modules obtained through the usual null model yield a Q value of 0.49 and are represented in Fig 3. Table 3 shows the different elements comprising each module. These modules are interpreted as the following morpho-functional complexes: (1) hip complex, tied with hip movements such as rotation, protraction, and flexion (green nodes in Fig 3); (2) thigh complex, including most of the nodes of the stylopod. This module is related to the distal  Table 2. The schematic representation of the involved anatomical elements is provided below. https://doi.org/10.1371/journal.pone.0177819.g003 The anatomical network of frog hindlimb movements, including adduction, retraction and protraction of the femur and the knee (black nodes in Fig 3); (3) shank complex, embracing the elements associated to the tibiofibula and fibulare, with associated muscles (blue nodes in Fig 3); (4) calf complex, including the prominent flexor digitorum communis and its origin fibrous knots (cian nodes in Fig 3); and (5) foot complex, embracing the fibulare and the multiple bones of the autopodium in addition to one extensor element and to the flexor surface of the foot (red nodes in Fig 3).

Discussion
We accept a basic statement of our work in that there is a correspondence between the topographical arrangements of elements and the network topology derived from their physical connections. This finding reinforces the notion that modeling the hindlimb as a network is a valuable resource for testing an old hypothesis from the field of anatomy (e.g. principe des  connexions about homologies [67]) and for studying scarcely explored ideas to date, such as tensegrity properties [68]. Bones are among the most central nodes, actually overpassing the fibrous knots. This result contradicts our first hypothesis, which considered the fibrous elements of connective tissue to be more prominent in this regard. Consequently, bones play a leading role on the integrity of the system. The fact that bones exhibit higher centrality values than elements such as fibrous knots is worth remarking, since the latter have been classically considered to present a main role connecting bones and muscles. Our network analysis, however, tells a different story in which bones are the frame connecting muscles and fibrous entities. A similar result was obtained by Diogo et al. [43] in relation to the human hindlimb. These results indicate that the main role of fibrous knots probably resides in their biomechanical properties, which would be much more important than their role as connectors. In fact, most of these fibrous knots tend to undergo metaplasic processes generating fibrocartilage, or even bones (e.g. the fibrocartilaginous patella in the knee knot [69], the plantar sesamoid in the plantar aponeurosis [70], and the cartilaginous sesamoid in the ligamentun calcanei [70]) as a response to the strong mechanical stress to which they are exposed [71]. This is another argument supporting their consideration as particular nodes in our network. Interestingly, in a recent work of Taglaferri et al. [72], the concept of bone-muscle unit is studied once again, without consideration of the tendinous system, thus reinforcing the notion of tendons as biological devices to handle the energetic intake needed to perform body movements, rather than as elements connecting other anatomical structures. Other authors [73,74] also concluded that one of the main functions of tendons is to facilitate joint motion, which exceeds the contractile capacities of the muscles.
There is a certain degree of correlation between osteological elements ordered by betweenness and the ontogenetic chronology of the hindlimbs in tetrapods. This is worth remarking, because it calls for a model of preferential attachment behind the organization of the network [75]. In other words, the osteological nodes which are first to come up are more likely to be connected with other elements that progressively arrive, and are thus prone to capture new connections, increasing their betweenness. In the ontogeny of the anuran hindlimb, the first osteological element to mature (chondrogenesis and osteogenesis) is the femur, followed by the tibiofibula and finally, by the autopodial elements [63]. More precisely, the elements with the higher level of betweenness, femur and tibiofibula, are the first to differentiate in the hindlimb. These elements have a pioneer role in determining a topographical differentiation inside the limb bud [63,76]. Interestingly, the acetabulum exhibits a lower level of betweenness than the long bones of the hindlimb, a rather counter intuitive result when considering the proximal-distal ontogenetic differentiation of the hindlimb. However, Manzano et al. [63] showed that the mesenchimal condensations of the hindlimb long bones originate in their diaphyses when the acetabulum is still in the interzone stage and when no joint cavity is visible, providing evidence of its late differentiation. Indeed, the level of betweenness of the pelvic girdle elements display the same phenomenon: they all exhibit lower level of betweenness than the hindlimb, in a consistent pattern with their ontogenetic sequence. Thus, our data reinforce the idea that the limb is incorporated as an almost complete framework for the development of the pelvic girdle. This topographical arrangement, together with the location of the tendon primordia as described by Manzano et al. [63], likely collaborates in the organization of the position of the muscular complexes [76].
The modularity analysis reveals a clustered organization of the tendon-musculo-skeletal network, supporting the implicit assertion of our second hypothesis. The modular arrangement adjusts to a topographical scheme of proximal-distal organization. Although this conclusion might seem obvious, it should be noted that other results, involving the mixing of distal and proximal elements, could have been obtained. Following this interpretation, the modular arrangement shows evidence of developmental and functional constraints in the structural organization of the frog hindlimb, sustaining thus our third hypothesis. Modules, interpreted as sets of densely connected elements, can be considered as an epiphenomenon comparable to that of anatomical modules from the evo-devo perspective. A clear example of this kind of reciprocal illumination includes the pervasive assignation of the autopodium to a module in its own right [77,78]. The modular construction of vertebrate appendages has been a major factor throughout their evolution. Shubin & Davis [79] conceive modules as existing at several hierarchical levels in limbs, from the whole organ to combinations of their endoskeletal constituents.
Two layers of nodes (skeletal and muscular) seem to exist, in which the elements of the latter are largely unconnected. If the single exception of muscle-to-muscle link is disregarded, the network as a whole would subsume into a type of graph known as quasi-bipartite graph. In a quasi-bipartite graph there are two sets of nodes: one partition with connections among them and another set of nodes with no connection. The adequacy of the quasi-bipartite model for musculoskeletal networks deserves a more in-depth consideration because of its potential impact for modularity studies. Given the explicit dependence of modularity upon the null model, it is clear that the specific choice of the null model largely influences modularity [80]. In fact, a quasi-bipartite graph calls for a different null model than the usual to study modularity, in which the constraints of the absence of a connection for the disjoint set of vertices should be considered. Therefore, the development of a more realistic null model for assessing the modular structure in musculoskeletal networks is of urgent need. For instance, P must assign zero likelihood to edges between vertices of the same independent set of vertices, precluding such edges in the null model. Thus, pairs of vertices from the naturally disjoint set should not receive a modularity penalty.
Jumping is a primary locomotion mode in frogs [47,81,82] and particularly in Leptodactylus latinasus, an outstanding jumping frog [83]. This remarkable and distinctive function can be disaggregated into functional sublevels, namely rotation, adduction, articulation, extension, etc. The coordinated action performed by modules within the hindlimb is thought to maximize jumping [84]. This is to be expected in the context of evolutionary biology, since modules are considered as autonomous anatomical parts, influenced by local factors and integrated through common developmental factors [40]. Functionally, there is a correlation between the proximal-to-distal succession of modules and the progressive recruitment of elements involved in joint motion during jumping. In a typical jump of a frog, the iliosacral and hip joint movements start almost simultaneously, followed by the knee joint and finally by the ankle joint [84,85]. This sequential movement of joints assures an extended acceleration of the center of mass [84].
The two elements that form the ilium -the iliac shaft and the acetabular portion-, belong to two different modules. One of them denotes the relationship between the pelvic girdle and the body, while the other exhibits a stronger relationship with the proximal elements of the thigh, i.e. the link between the pelvic girdle and the hindlimb. This separation of the same bone in two modules could be interpreted from a functional perspective, since the elongated ilial shafts of anurans are specifically modified for jumping [86]; the acetabular portion, instead, acts as an attachment element with the hindlimbs, through the joint with the head of the femur, in a similar way than in other tetrapods.
The thigh module is mostly related to the flexion and extension of the shank at the knee [87] and is marked by a relatively high percentage of fibrous elements (fascia latae, knee knot, tenuissimus it knot). In the network layout, guided by topological information alone, this module seems to bisect the proximal and distal pools of anatomical elements.
The flexor digitorum communis muscle is clearly the strongest muscle of the hindlimb, and is the main muscle used by frogs in jumping [88]. It is the only muscular element which comprises the calf module, along with its fibrous origin knots. The outer links of this module correspond to the distal part of the femur and the aponeurosis plantaris. This last edge is the Achilles tendon, responsible of the storage of elastic energy for the amplification of power output during jumping [89]. Further studies are needed to assess whether these modules are also recovered in the network of hindlimbs of other species of jumping tetrapods.
The anuran hindlimb presents an additional functional segment composed by the tibiale and fibulare [90,91], which is also incorporated in musculoskeletical modeling [92]. These bones are usually short bones belonging to the autopodium, but in anurans they elongate to form the additional functional segment mentioned above. In our analysis, the tibiale is included in a common module, also integrated by the tibiofibula. Therefore, at least in relation to this bone, our network reinforces the notion of zeugopodization experienced by certain mesopodial bones that receive signals coming from the zeugopodial area to start elongation [93].
The most distal module comprises the many bones of the foot, the aponeurosis plantaris and two muscles which are responsible of the extension of the digits and of the dorsiflexion and pronation of the foot [88]. Diogo et al. [43] included more elements in the autopodial segment than those included in our analysis, and they were able to outline six modules in their study. In our analysis, the elements of the foot are grouped in a single module. Our data show a lower degree of resolution and the discrepancy in the number of modules might be a matter of scale. To conclude, modular architecture seems to be a successful organization that provides the building blocks over which evolution forges the many different functional specializations that organisms can exploit.