Growing networks-backup.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:f87e8922f6540596a32443d0b7842c1db89b28841872f95e1345250e99ca2b1e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# A growing brain network model\n",
      "*Roberto Toro, Septembre 2014*"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It has been reported that the topology of brain networks optimises wiring and maximise efficiency (Bullmore and Sporns, 2012). This optimisation has been regarded as the fruit of evolution, a quality that may have been acquired through millions of years of adaptive natural selection. The null hypotheses against which optimisation is judged -- random networks or regular lattices -- are, however, rather irrealistic. Brain networks are the result of natural developmental phenomena, the growth of dendritic trees and axons in a substrate that is itself undergoing massive growing. Newly generated neurones have been observed to send axons to subcortical structures since the earliest stages, well before reaching the cortical plate. Some of these connection may be pruned at later stages, but some of them will form the basic pattern of brain connections that will be still present in the adult. Early connections will be formed when the brain volume is just a few cm3, but will mature into the connections observed in the adult, when the brain volume will be more than 1,000 times larger. Producing a random graph topology or a regular lattice topology would require developmental processess incredibly complex."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here we propose a simple model of brain network generation intended to serve as a more realistic null hypothesis. The model is based on just two very basic hypotheses: (1) the number of neurones increases through time, (2) the probability of two neurones forming connections is inversely proportional to their distance. As the brain model grows, short-distance connections formed among regions that were initially close, become long-distance connections. There will then be a link between the time of formation of a connection and its length: earlier connections will on average longer than the later ones.\n",
      "\n",
      "We start by analysing the topological properties of these graphs, their small-worldness, modularity, average short path length, their efficiency, and then we use them as null hypothesis to evaluate real structural and functional connectivity graphs from human neuroimaging data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import networkx as nx\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "G=nx.barabasi_albert_graph(100,5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pos=nx.graphviz_layout(G,prog=\"neato\")\n",
      "# nx.draw_random(G)\n",
      "nx.draw(G,pos)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"avg. clustering:\",np.mean(nx.clustering(G).values())\n",
      "print \"avg. deg. centrality:\",np.mean(nx.degree_centrality(G).values())\n",
      "print \"diameter:\",nx.diameter(G)\n",
      "print \"avg. shortest path length:\",nx.average_shortest_path_length(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "avg. clustering: 0.195603255745\n",
        "avg. deg. centrality: 0.0959595959596\n",
        "diameter: 4\n",
        "avg. shortest path length: 2.21292929293\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Brain networks do not show significantly different topological properties than those generated by our model. In population genetics, the two main forces driving evolution are selection and natural drift. Selection captures the adapatative modification of the phenotype, whereas drift captures the random effects due to environmental constraints such as distance, nanani, nanana. The current phenotype of a species is not only the result of adaptive selection, but more exactly the phenotype that has slip through the nets of adversity. Spandrels, Gould & Lewontin describe this nicely, use what you have at hand. Modern human brains have existed for some 130k years, however, the first traces of what we consider human culture, art, pottery, are much more recent, maybe less than 60k years. It is probable that for some 70k years, human had a modern brain topology but nature had yet to discover what it was good for."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import print_function, division\n",
      "import networkx as nx\n",
      "\n",
      "\n",
      "def global_efficiency(G, weight=None):\n",
      "    \"\"\"Return the global efficiency of the graph G\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    G : NetworkX graph\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    global_efficiency : float\n",
      "\n",
      "    Notes\n",
      "    -----\n",
      "    The published definition includes a scale factor based on a completely\n",
      "    connected graph. In the case of an unweighted network, the scaling factor\n",
      "    is 1 and can be ignored. In the case of a weighted graph, calculating the\n",
      "    scaling factor requires somehow knowing the weights of the edges required\n",
      "    to make a completely connected graph. Since that knowlege may not exist,\n",
      "    the scaling factor is not included. If that knowlege exists, construct the\n",
      "    corresponding weighted graph and calculate its global_efficiency to scale\n",
      "    the weighted graph.\n",
      "\n",
      "    Distance between nodes is calculated as the sum of weights. If the graph is\n",
      "    defined such that a higher weight represents a stronger connection,\n",
      "    distance should be represented by 1/weight. In this case, use the invert_\n",
      "    weights function to generate a graph where the weights are set to 1/weight\n",
      "    and then calculate efficiency\n",
      "\n",
      "    References\n",
      "    ----------\n",
      "    .. [1] Latora, V., and Marchiori, M. (2001). Efficient behavior of\n",
      "       small-world networks. Physical Review Letters 87.\n",
      "    .. [2] Latora, V., and Marchiori, M. (2003). Economic small-world behavior\n",
      "       in weighted networks. Eur Phys J B 32, 249-263.\n",
      "\n",
      "    \"\"\"\n",
      "    N = len(G)\n",
      "    if N < 2:\n",
      "        return 0    # facilitates calculation of local_efficiency although\n",
      "                    # could reasonably raise nx.NetworkXUnfeasible or\n",
      "                    # nx.NetworkXPointlessConcept error instead and force\n",
      "                    # testing to occur in local_efficiency\n",
      "\n",
      "    inv_lengths = []\n",
      "    for node in G:\n",
      "        if weight is None:\n",
      "            lengths = nx.single_source_shortest_path_length(G, node)\n",
      "        else:\n",
      "            lengths = nx.single_source_dijkstra_path_length(G, node,\n",
      "                                                            weight=weight)\n",
      "\n",
      "        inv = [1/x for x in lengths.values() if x is not 0]\n",
      "        inv_lengths.extend(inv)\n",
      "\n",
      "    return sum(inv_lengths)/(N*(N-1))\n",
      "\n",
      "\n",
      "def local_efficiency(G, weight=None):\n",
      "    \"\"\"Return the local efficiency of each node in the graph G\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    G : NetworkX graph\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    local_efficiency : dict\n",
      "       the keys of the dict are the nodes in the graph G and the corresponding\n",
      "       values are local efficiencies of each node\n",
      "\n",
      "    Notes\n",
      "    -----\n",
      "    The published definition includes a scale factor based on a completely\n",
      "    connected graph. In the case of an unweighted network, the scaling factor\n",
      "    is 1 and can be ignored. In the case of a weighted graph, calculating the\n",
      "    scaling factor requires somehow knowing the weights of the edges required\n",
      "    to make a completely connected graph. Since that knowlege may not exist,\n",
      "    the scaling factor is not included. If that knowlege exists, construct the\n",
      "    corresponding weighted graph and calculate its local_efficiency to scale\n",
      "    the weighted graph.\n",
      "\n",
      "    References\n",
      "    ----------\n",
      "    .. [1] Latora, V., and Marchiori, M. (2001). Efficient behavior of\n",
      "       small-world networks. Physical Review Letters 87.\n",
      "    .. [2] Latora, V., and Marchiori, M. (2003). Economic small-world behavior\n",
      "       in weighted networks. Eur Phys J B 32, 249-263.\n",
      "\n",
      "    \"\"\"\n",
      "    if G.is_directed():\n",
      "        new_graph = nx.DiGraph\n",
      "    else:\n",
      "        new_graph = nx.Graph\n",
      "\n",
      "    efficiencies = dict()\n",
      "    for node in G:\n",
      "        temp_G = new_graph()\n",
      "        temp_G.add_nodes_from(G.neighbors(node))\n",
      "        for neighbor in G.neighbors(node):\n",
      "            for (n1, n2) in G.edges(neighbor):\n",
      "                if (n1 in temp_G) and (n2 in temp_G):\n",
      "                    temp_G.add_edge(n1, n2)\n",
      "\n",
      "        if weight is not None:\n",
      "            for (n1, n2) in temp_G.edges():\n",
      "                temp_G[n1][n2][weight] = G[n1][n2][weight]\n",
      "\n",
      "        efficiencies[node] = global_efficiency(temp_G, weight)\n",
      "\n",
      "    return efficiencies\n",
      "\n",
      "\n",
      "def average_local_efficiency(G, weight=None):\n",
      "    \"\"\"Return the average local efficiency of all of the nodes in the graph G\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    G : NetworkX graph\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    average_local_efficiency : float\n",
      "\n",
      "    Notes\n",
      "    -----\n",
      "    The published definition includes a scale factor based on a completely\n",
      "    connected graph. In the case of an unweighted network, the scaling factor\n",
      "    is 1 and can be ignored. In the case of a weighted graph, calculating the\n",
      "    scaling factor requires somehow knowing the weights of the edges required\n",
      "    to make a completely connected graph. Since that knowlege may not exist,\n",
      "    the scaling factor is not included. If that knowlege existed, a revised\n",
      "    version of this function would be required.\n",
      "\n",
      "    References\n",
      "    ----------\n",
      "    .. [1] Latora, V., and Marchiori, M. (2001). Efficient behavior of\n",
      "       small-world networks. Physical Review Letters 87.\n",
      "    .. [2] Latora, V., and Marchiori, M. (2003). Economic small-world behavior\n",
      "       in weighted networks. Eur Phys J B 32, 249-263.\n",
      "\n",
      "    \"\"\"\n",
      "    eff = local_efficiency(G, weight)\n",
      "    total = sum(eff.values())\n",
      "    N = len(eff)\n",
      "    return total/N\n",
      "\n",
      "\n",
      "def invert_weights(G, weight='weight'):\n",
      "    \"\"\"Return a graph where the weight of each edge is 1 over the weight of\n",
      "       the corresponding edge in graph G\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    G : NetworkX graph\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    Ginv : NetworkX graph\n",
      "\n",
      "    Notes\n",
      "    -----\n",
      "    This function is meant to address the case where weights represent the\n",
      "    \"strength\" of the connection rather than the distance between nodes. In\n",
      "    this case, the distance would be considered to be 1 over the weight so\n",
      "    that \"weak\" connections have correspondingly \"long\" distances between\n",
      "    nodes.\n",
      "\n",
      "    \"\"\"\n",
      "    Ginv = nx.create_empty_copy(G)\n",
      "    for (n1, n2) in G.edges():\n",
      "        dist = 1/G[n1][n2][weight]\n",
      "        Ginv.add_edge(n1, n2, {weight: dist})\n",
      "\n",
      "    return Ginv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"globale efficiency:\",global_efficiency(G))\n",
      "print(\"average local efficiency:\",average_local_efficiency(G))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "globale efficiency: 0.496666666667\n",
        "average local efficiency: 0.328854003039\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "le=local_efficiency(G)\n",
      "nx.draw(G,pos,node_color=[value for key,value in le.iteritems()],cmap = plt.cm.Blues)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import mayavi.mlab as mlab\n",
      "pos=nx.spring_layout(G,dim=3)\n",
      "xyz=np.array([pos[v] for v in sorted(G)])\n",
      "colors=[value for key,value in le.iteritems()]\n",
      "\n",
      "mlab.figure(1, bgcolor=(0, 0, 0))\n",
      "mlab.clf()\n",
      "pts = mlab.points3d(xyz[:,0], xyz[:,1], xyz[:,2],\n",
      "                    colors,\n",
      "                    scale_factor=0.05,\n",
      "                    scale_mode='none',\n",
      "                    colormap='Blues',\n",
      "                    resolution=20)\n",
      "\n",
      "pts.mlab_source.dataset.lines = np.array(G.edges())\n",
      "tube = mlab.pipeline.tube(pts, tube_radius=0.001)\n",
      "mlab.pipeline.surface(tube, color=(0.8, 0.8, 0.8))\n",
      "#mlab.savefig('mayavi2_spring.png')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 34,
       "text": [
        "<mayavi.modules.surface.Surface at 0x1100d5e90>"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}