[Pybrainsim-activity] SF.net SVN: pybrainsim:[127] trunk/src/PyBrainSim.py

[<rg...@us...>]

Revision: 127
          http://pybrainsim.svn.sourceforge.net/pybrainsim/?rev=127&view=rev
Author:   rgoj
Date:     2009-09-13 19:24:15 +0000 (Sun, 13 Sep 2009)

Log Message:
-----------
* Added an example cortical area simulation with a neural mass model from Jansen and Rit 1995.

This currently isn't programmed with classes, so it will need to be rewritten to allow neural mass model specification, but it shows that it is easy to use such models and is a good starting point for writing the classes.

Modified Paths:
--------------
    trunk/src/PyBrainSim.py

Modified: trunk/src/PyBrainSim.py
===================================================================
--- trunk/src/PyBrainSim.py	2009-09-11 10:30:28 UTC (rev 126)
+++ trunk/src/PyBrainSim.py	2009-09-13 19:24:15 UTC (rev 127)
@@ -178,6 +178,88 @@
     exampleHead = Head()
     exampleHead.displayHead()
 elif userChoice == 8:
-    print("Not implemented yet...")
+    # This class will produce noise that will be fed into the modelled
+    # cortical area as input.
+    class NeuralNoise:
+        def __init__(self, timeSpan, dtime):
+            self.noise = []
+            self.dtime = dtime
+            self.timePoints = range( int(timeSpan//dtime+1) )
+            # Random noise distributed evenly between 120 and 320
+            for i in self.timePoints:
+                self.noise.append( 120+200*random() )
+        def getNoise(self, time):
+            # Searching for the appropriate time point
+            for i in self.timePoints:
+                if time < self.timePoints[i]*self.dtime:
+                    # Interpolating between the randomly generated values
+                    return self.noise[i] + \
+                           (self.noise[i+1]-self.noise[i]) \
+                           * ((self.timePoints[i]*self.dtime-time)/self.dtime)
+    
+    # The sigmoid function
+    def sigmoidFunction(v):
+        """Constants taken from Jansen and Rit 1995, p. 360"""
+        cEzero = 2.5
+        cR = 0.56
+        cVzero = 6
+    
+        """Sigmoid function taken from Jansen and Rit 1995, p. 360"""
+        return 2*cEzero / ( 1 + numpy.exp( cR*( cVzero - v ) ) )
+    
+    # Defining the equations
+    def f(t, y, noiseObject, cC):
+        cA=3.25
+        ca=100.0
+        cB=22.0
+        cb=50.0
+        cC1=cC
+        cC2=0.8*cC1
+        cC3=0.25*cC1
+        cC4=0.25*cC1
+        noise = noiseObject.getNoise(t)
+        #print("Time: " + str(t) + "   Noise: " + str(noise))
+        return [y[1], cA*ca*( sigmoidFunction(y[2] - y[4]) ) - 2*ca*y[1] -
+                ca**2*y[0], y[3], cA*ca*( noise + cC2*sigmoidFunction(cC1*y[0]) ) -
+                2*ca*y[3] - ca**2*y[2], y[5], cB*cb*( cC4*sigmoidFunction(cC3*y[0])
+                                                ) - 2*cb*y[5] - cb**2*y[4]]
+    
+    # How many seconds should be modelled and how accurately
+    timeSpan = 3
+    dtime = 0.001
+    firstPointToDisplay = 1000
+
+    # A noise object that will be used as input into the equations
+    someNoise = NeuralNoise(timeSpan+1, dtime*10)
+    
+    # Initial conditions
+    y0, t0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 0
+
+    # Preparing for simulations
+    r = []
+    recording = [ [], [],  [],  [],  [],  [] ]
+    parameter = [ 68, 128, 135, 270, 675, 1350 ]
+    
+    # Differential equation integration
+    for i in range(len(parameter)):
+        # Preparing the integration
+        r.append(ode(f).set_integrator('vode', method='bdf'))
+        r[i].set_initial_value(y0, t0)
+        r[i].set_f_params( someNoise, parameter[i] )
+        # Integrating
+        print("Performing simulation " + str(i+1) + "/" + str(len(parameter)))
+        while r[i].successful() and r[i].t < timeSpan:
+            r[i].integrate(r[i].t+dtime)
+            recording[i].append(r[i].y[2] - r[i].y[4])
+
+    # Show results of simulations
+    from pylab import *
+    for i in range(len(parameter)):
+        subplot( len(parameter), 1, i+1)
+        if i==0:
+            title("Simulations based on Fig. 3 from Jansen and Rit 1995")
+        ylabel( "C = " + str(parameter[i]) )
+        plot(recording[i][firstPointToDisplay:])
+    show()
 else:
     print("No such option unfortunately...")


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.