Mathematica Package
Mathematica Package Listing
1 (* Package: PhaseFunction *)
2 (* Author: Ramana Dodla, 28 August 2009 *)
3 (* Reference: Ramana Dodla and Charles J. Wilson. Biophysical Journal 98:L05-L07, 2010 *)
4 (* Comment: User should implement methods, if (s)he wishes to avoid errors for impossible lag (tau) values *)
5
6 BeginPackage["PhaseFunction`"]
7 PhaseFunction::usage="This package provides these three functions:
8 1. Psi[A_,B_,tau_] that finds the phase function between spike trains A and B at a given time lag tau,
9 2. PhaseDispC[A_,B_] that finds the list c_i, and
10 3. PhaseDispPairs[A_,B_] that finds (alpha_i,beta_i)
11 Example:
12 A={0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024}
13 ListPlot[Table[{i*0.001, Psi[A,A,i*0.001]}, {i, 0, 400}]]
14 "
15 Psi::usage="Psi[spiketrain1, spikerain2, tau]"
16 PhaseDispC::usage="PhaseDispC[spiketrin1, spiketra2] returns the list c_i, the drift parameter"
17 PhaseDispPairs::usage="PhaseDispPairs[spiketrain1, spiketain2] returns (alpha_i, beta_i) pairs"
18 Begin["`Private`"]
19
20 (* Compute phase function, psi *)
21 Psi[A_, B_, tau_] := Module[
22 {c = {}, norm, thetadrift},
23 c = PhaseDispC[A, B + tau];
24 thetadrift = Mean[Abs[c[[All]]]];
25 norm = Sqrt[2.0];
26 thetadrift / norm
27 ];
28 PhaseDispC[A_, B_] := Module[
29 {mc, md, f1, f2, c = {}, Pairs = {}},
30 Pairs = PhaseDispPairs[A, B];
31 mc = Mean[Pairs[[All, 1]]];
32 md = Mean[Pairs[[All, 2]]];
33 f1 = mc * mc + md * md;
34 f2 = Sqrt[f1];
35 c = (f1 - Pairs[[All,1]] * mc
36 - Pairs[[All,2]] * md) / f2
37 ];
38 PhaseDispPairs[A_, B_] := Module[
39 {i, j, g, d, T, Pairs = {}},
40 i = 2;
41 j = 2;
42 T = Max[{A[[1]], B[[1]]}];
43 While[ i <= Length[A] && j <= Length[B],
44 If[ A[[i]] >= B[[j]],
45 If[A[[i-1]] >= B[[j]],
46 T = A[[i-1]];
47 j = j+1,
48 g = (B[[j]] - T) / (A[[i]] - A[[i-1]]);
49 d = (B[[j]] - T) / (B[[j]] - B[[j-1]]);
50 AppendTo[Pairs, {g, d}];
51 T = B[[j]];
52 j = j + 1
53 ],
54 If[B[[j-1]] >= A[[i]],
55 T = B[[j-1]];
56 i = i + 1,
57 g = (A[[i]] - T) / (A[[i]] - A[[i-1]]);
58 d = (A[[i]] - T) / (B[[j]] - B[[j-1]]);
59 AppendTo[Pairs, {g, d}];
60 T = A[[i]];
61 i = i + 1
62 ];
63 ];
64 ];
65 Pairs
66 ];
67
68 End[]
69 EndPackage[]