WO2002003327A2

WO2002003327A2 - Method and system for analyzing multi-dimensional data

Info

Publication number: WO2002003327A2
Application number: PCT/US2001/021032
Authority: WO
Inventors: Rafael Yuste; Vikram S. Kumar; Robert C. Froemke
Original assignee: The Trustees Of Columbia University In The City Of New York
Priority date: 2000-06-29
Filing date: 2001-06-29
Publication date: 2002-01-10
Also published as: AU2001271761A1; US20040015310A1; WO2002003327A3; US6823286B2

Abstract

Disclosed is a method for analyzing a sequence of data arrays. A selection of at least one type of region of interest and at least one region of interest from said data arrays is made. The sequence of data arrays are then transformed into a simplified data array. Events of interest in the selected regions of interest are then detected and stored in a second simplified data array. The data is then analyzed to determine relationships between the detected events of interest.

Description

METHOD AND SYSTEM FOR ANALYZING MULTI-DIMENSIONAL DATA

SPECIFICATION

RELATED APPLICATION This application claims priority from U.S. provisional application No. 60/214,914 filed on June 29, 2000, which is incorporated by reference herein in its entirety

BACKGROUND OF INVENTION The present invention relates to analyzing and inteφreting multidimensional datasets. Examples of such datasets include optical recordings of neuronal cell slice fluorescenceand differences in expression levels of multiple genes within a population of patients or subjects.

It is often desirable to understand the relationship of various events occurring within such a multidimensional dataset. For example, various neurons in a neuronal cell slice may exhibit spontaneous activity in a time series of optical images. It would be desirable to determine which, if any, group of neurons were ever coactive (i.e. active at the same time or at specifiec different times), were regularly coactive (i.e. coactive at multiple times over the period of observation), and which neuron, if any, consistently activates before or after another neuron activates. It would also be advantageous to know the statistical significance of the relationships between the various events. In other words, whether the correlation among the various events is stronger than would be expected from random activity.

SUMMARY OF THE INVENTION These and other advantages are achieved by the present invention which provides a method and system for analyzing a multi-dimensional dataset and for detecting relationships between various events reflected in the dataset. In an exemplary embodiment, a method is presented for analyzing a sequence of data arrays including selecting at least one type of region of interest and at least one region of interest for each type of region of interest chosen from said data arrays, and transforming the sequence of data arrays into a simplified data array with a first dimension equal to the number of selected regions of interest and a second dimension equal to the number of data arrays in the original sequence of data arrays. The simplified data array is then examined to detect events of interest in the regions of interest, and those events of interest are stored in a second simplified data array having the same dimensions as the first simplified data array, but the data in each element of the array is binary. The second simplified array is then analyzed to determine relationships between the events of interest and correspondingly, the regions of interest.

In one exemplary embodiment, analyzing includes plotting a portion or all of the data in the first simplified array to allow visual examination of the relationships between the activities of interest in various regions of interest. In another exemplary embodiment, the analysis step involves detecting events of interest that are coactive and determining whether the number of coactive events is statistically significant. This embodiment may include detecting all such coactive events (i.e. events where at least two regions of interest are active simultaneously), detecting instances where many regions of interest are coactive simultaneously, or detecting instances where two or more regions of interest are each active in a certain temporal relationship with respect to ont another (also referred to as coactivity).

In a further exemplary embodiment, the data analysis involves calculating a correlation coefficient between two regions of interest based on how often the regions of interest are coactive relative to how often the first region is active. A map of all such regions is displayed with lines between the regions having a thickness proportional to the correlation coefficient between the two regions.

Another exemplary embodiment includes plotting a cross-correlogram or histogram of events of interest in a particular region of interest with respect to events of interest in another region of interest, so that the histogram will reveal the number of times an event of interest in the first region of interest occurs a certain number of locations away from an event of interest in the second region of interest in the second simplified data array. The cross-correlogram can be plotted with respect to one region of interest, thus showing how many times an event of interest occurs before or after the occurance of another event of interest in the same region of interest.

Other exemplary embodiments include performing Hidden Markov

Modeling on the second simplified data array to determine a hidden Markov state sequence and displaying a cross-correlogram between events of interest occurring in one region of interest while that region is in one of the detected Markov states and performing a singular value decomposition on the first simplified data array.

BRTEF DISCRIPTION OF THE DRAWINGS For a more complete understanding of the present invention, reference is made to the following detailed description of a exemplary embodiments with reference to the accompanying drawings in which: Fig. 1 illustrates a flow diagram of a method in accordance with the present invention;

Fig. 2 illustrates an example sequence of data arrays for analyzing in accordance with the present invention;

Fig. 3 illustrates an example of a simplified data structure generated by the method of Fig. 1;

Fig. 4 illustrates a flow diagram of a method of analyzing data useful in the method of Fig. 1;

Fig. 5 illustrates a visual plot generated in accorance with the method of Fig. 1; Fig. 6 illustrates a cross-correlogram generated in accordance with the method of Fig. 1; and

Fig. 7 illustraes a correlation map generated in accordance with the method of Fig. 1. DESCRIPTION OF THE PREFERRED EMBODIMENT Referring to Fig. 1, there is shown a flow diagram represetning an exemplary method for analyzing a sequence of data arrays in accordance with the present invention. For purposes of this description, the sequence of data arrays operated upon corresponds to a series of two dimensional optical images of neuronal cell slices, such as a slice of brain tissue, captured at a fixed interval of time. Thus, the input data is a series of two dimensional data arrays, with each dimension corresponding to a spatial dimension and each element of the data array corresponding to fluorescence level or other optical measure of activity in the imaged neuronal cell slice. Each data array in the series corresponds to image data at a different instant in time, with the images taken at fixed intervals of time. The format of this input data will be discussed further herein with reference to Fig. 2. It will be understood that the present invention is not limited to such data. For example, the input data could correspond to expression levels of genes within a population of subjects. Other potential input data sets will be apparent to one of ordinary skill in the art.

In the exemplary embodiment, performance of the method is assisted by a general purpose computer adapted to operate the MAC-OS operating system and to interpret program code written in Interactive Data Language ("IDL") version 5.1 or later, developed by Research Systems, Inc. The IDL program code of the exemplary embodiment is appended hereto as Appendices A, B and C described further herein. Other operating systems and programming languages could be used to perform the steps of the exemplary embodiment without departing from the scope of the invention, and the modifications necessary to make such a change will be apparent to one of ordinary skill in the art.

In step 101, a selection of at least one type of region of interest of the input data arrays is made. This selection may be made manually, automatically based on certain criteria relevant to the particular data being analyzed, or predetermined. In the exemplary embodiment, there is only one relevant type of region of interest, and that is predetermined to be a neuron in the image data. If other data were to be analyzed, the relevant type of region of interest may be different and there may be multiple such types of regions.

In step 103, a selection of one or more regions of interest for each of the selected types of regions of interest is made. To further understand this step in the exemplary embodiment, reference is made to Fig. 2 where an example of an input data array of a sequence of two dimensional image data is shown. Data array 201 is a array of picture element (pixel) data corresponding to an image acquired at time t=l. As can be seen, in this data array, there are sixteen pixels 207, each representing the optical intensity at a corresponding spatial region of the original image. Although only sixteen pixels are shown here, it will be understood that actual data will have many more pixels, corresponding to the resolution of the imaging device. In the preferred embodiment, the imaging device is a charge coupled device (CCD) camera which is capable of converting an optical image into digital pixel data.

A fixed interval of time after the image array 201 is acquired, second image array 203 is acquired. The process of acquiring data images continues at fixed intervals until final image array 205 is acquired at time t=N. In this example, time t=l corresponds to the time the first image was take and time t=N corresponds to the time the last image was taken. The fixed time interval and number of images acquired (N) can be adjusted to maximize time resolution by reducing the fixed time interval. By reducing the time interval and increasing N, the amount of data collected will increase accordingly, which may lead to longer analysis times.

Returning to Fig. 1, the step of selecting a region of interest 103 consists of selecting an area of the image data where a neuron is present. For example, in Fig. 2, the circles 209 and 223 represent neurons in the original image data that are the regions of interest. Accordingly, the selecting step 103 consists of selecting the image data at pixel locations 211 and 217 as regions of interest. In the preferred embodiment, the step is performed manually by a user viewing an image corresponding to image data frame 201 or other image data frame and selecting the area where the neuron is visible. Because each image pixel data array in the input data sequence corresponds to the same spatial area over time, and because the neurons do not move over the time of observation, the pixels of interest in all subsequent images arrays can also be identified. Thus, in Fig. 2, it can be seen that pixel 219 in image pixel array 203 and pixel 221 in image pixel array 205 represent the same neuron as pixel 217 in image array 201. Similarly, pixels 211, 213, and 215 all represent the same neuron of interest. In step 105, the data corresponding to the selected regions of interest are tranformed into a simplified data array. This simplified data array contains only the pixel intensity data from the original sequence of pixel data arrays that corresponds to the regions of interest. An example of such a simplified data array is shown in Fig. 3. There, each data row corresponds to one of the selected regions of interest from the original sequence of pixel data arrays, while each column represents a different image pixel array in the original sequence of image pixel arrays. Each element in the simplified data array contains the intensity value of pixels in the relevant region of interest at a particular time. Thus, in step 105, the intensity data for pixel 211 in image pixel array 201 is placed in the portion of simplified data array 300 designated as array element 301, while pixel intensity data 217 is placed in array element 307. Similarly, pixel intensity data 213 and 219 are placed in array elements 303 and 309 respectively. The process continues until pixel intensity data 215 and 221 are placed in elements 305 and 311 respectively. It will be understood that if more than two neurons were selected in step 103, simplified array 300 would have more than 2 rows. If a neuron selected in step 103 spans more than one pixel position in the original sequence of image arrays, the average pixel intensity of the pixels corresponding to each selected neuron can be stored in the elements of simplified data array 300. If all of the data in the original sequence of data arrays corresponds to regions of interest, the transformation step would consists only of storing the input data into a data array having the dimensions described and no data from the original sequence need be ignored or discarded during the transformation step.

In step 107, events of interest in the simplified data array 300 are detected. In the exemplarly embodiment an event of interest is detected by calculating a statistical mean and standard deviation for all pixel intensity data corresponding to a particular region of interest. Thus, where the pixel intensity data is contained in the simplified array 300, a mean and standard deviation is calculated for all data pixel intensity in each row of the simplified array. An event is then detected where the pixel intensity data for a particular region of interest exceeds the mean for all data in the region of interest by a predetermined number of standard deviations. If activity were defined by a drop in intensity rather than an increase in intensity, the event could be detected by examining the pixel intensity data in a certain region of interest for an entry where the intensity is less than the mean for all data in the region of interest by a predeterminied number of standard deviations. The number of standard deviations may be entered by a user before the calculations are preformed, or a default number may be used, such as two or three. In this fashion, the method will detect those instances in time where the pixel intensity is much higher than the average intensity, thus suggesting neuron activity has occurred.

In another exemplary embodiment, an event is detected by looking for pixel intensities that exceed previous pixel intensities of the same regions of interest by a threshold amount. Thus, for example, if the pixel intensity stored in data element 309 exceeded the pixel intensity stored in data element 307 by a certain percentage, an event is said to have occurred at the time corresponding to pixel position 307 (i.e. t=l). Again, if an event were indicated by a drop in intensity rather than an increase, the detection step would involve looking for pixel intensity that are less than previous pixel intensities of the same region of interest by the threshold amount. The threshold amount can be specified by a user before the calculations are performed, or a default number can be used such as twenty percent. The detection can occur over many time periods, for example, the data corresponding to the image taken at time t=6 could be compared to the image taken at time t=l to see if an increase beyond the threshold amount has occurred. This would be useful to detect events that occur gradually over time rather than relatively instantaneously.

In step 109, the results of detection step 107 are stored in a second simplified array. For this purpose, the second simplified array is identical to the first simplified array illustrated in Fig. 3; however, the data stored in the second simplified array is binary rather than pixel intensity values. Thus, the entries in the second simplified array would be 1 or 0 (or yes or no), corresponding to whether an event of interest occurred in that region of interest at the corresponding time.

In step 111, the stored data is analyzed. In the preferred embodiment, the data is analyzed to determine whether various neurons are correlated (i.e. whether they are coactive), the strength of those correlations (i.e. how often they are coactive relative to how many times each neuron or one of the neurons is active), how significant the correlations are (i.e. whether the correlation is stronger than would be expected if from a random data set) and the behavior of the entire neuron population. In the exemplary embodiment, the data is analyzed by plotting at least a portion of the data contained in the first simplified data array 300. For example, pixel intensity for one neuron can be plotted over time. Pixel intensities for all neurons could also be plotted over time, either in separate plot windows or superimposed on the same plot window. Additionally, the pixel intensities for all neurons could be averaged and plotted over time to show global behavior of the systems. Fig. 5 illustrates three possible plots 501, 502, 503 of pixel intensity over time.

In another exemplary embodiment illustrated in Fig. 4, the data in the second simplified array is analyzed to determine the number of coactive events in the dataset and the statistical significance of those events. In step 401, a random distribution of neuron activity is generated. The random data is generated by shifting the data in each row of the second simplified array by a random amount. In step 403, the number of coactive events in the random dataset is counted. This process is repeated numerous time to generate a random distribution. The number of random trials may be set by the user or a default number of random trials may be conducted, such as 1000.

Counting coactive events for this purpose means counting all instances where two neurons are coactive. Coactive events for this purpose means events of interest that occurred in two neurons at the same time, or within a specified number of time intervals from each other. Thus, if the specified number of time intervals is one, then if a event occurred in neuron 1 at time t=l 301 and an event occurred in neuron 2 at time t=2 309, those events would be considered coactive. The time interval may be specified by a user before coactive events are counted, or may be a default setting such as two time intervals.

Once the random trials have been completed and a random distribution of coactive events generated, the actual number of coactive events in the data is calculated in step 405 using the same counting methodology was used to count coactive events in the random trials. The actual number of coactive events is then superimposed on a plot of the random distribution. The statistical significance of the coactive events is determined in step 407 by calculating the area under the distribution curve to the right of the number of actual coactive events in the data. This result, termed the "p-value" represents the probability that the number of detected coactive events in the actual data is produced by a random neuron activity.

In a further exemplary embodiment, a random distribution of activity is generated as previously described, except the only coactive events that are counted in steps 403 and 405 are those where a predetermined number of neurons are coactive. The predetermined amount of coactive neurons may be specified by a user or a predetermined default value such as four may be used. Additionally, it may be specified whether exactly that many coactive events must be present or at least that many coactive events must be present to be considered a coactive event for counting. Thus, the embodiment allows instances of multiple neurons active simultaneously (rather than simply two neurons active simultaneously) to be counted and the statistical significance of that number to be reported. In this exemplary embodiment, the random distribution and actual number of coactive events are plotted. The statistical significance of the actual number of coactive events is calculated using the formula: C_rand N_rand where C_rand is the number of random trials that resulted in more coactive matches than the actual data set and N_rand is the total number of random trials used to generate the random distribution, and is reported to a user. Additionally, a network map may be plotted showing the spatial locations of the coactive neurons with lines between those that were coactive. The spatial locations of the active neurons were those specified by the user during neuron selection step 103 in Fig. 1. In a still further exemplary embodiment, a random distribution of neuron activity is generated as previously described except the only coactive events that are counted in steps 403 and 405 are those where at least two neurons are active a predetermined number times throughout the dataset. The number of times the two or more neurons must be active can be specified by a user or a default number such as two may be used. In this exemplary embodiment, the random distribution and actual number of coactive events are plotted. The statistical significance of the actual number of coactive events is calculated using the formula: C_rand/N_rand where C_rarιd is the number of random trials that resulted in more coactive matches than the actual data set and N_rand is the total number of random trials used to generate the random distribution, and is reported to a user. Additionally, a network map may be plotted showing the spatial locations of the coactive neurons with lines between those neurons that were coactive the specified number of times. The spatial locations of the active neurons were those specified by the user during neuron selection step 103 in Fig. 1.

In yet another preferred embodiment, a correlation map is plotted. To plot the correlation map, a correlation coefficient array is first generated for all of the neurons. The correlation coefficients are defined as C(A,B)=number of times neuron A and B are coactive divided by the number of times neuron A is active. For this purpose, coactive means active at the same time, or within a specified number of time intervals of each other. The number of time intervals may be specified by a user or a default number such as one time increment may be used. The number of correlation coefficients will be equal to the square of the number of neurons selected in step 103 in Fig. 1. A correlation map is then drawn consisting of a map of all active neurons with lines between each pair of neurons having a line thickness proportional to the correlation coefficient of those two neurons. An example of such a correlation map is illustrated in Fig. 7. There, the thiclcness of line 707 is proportional to the magnitude of the correlation coefficent for neurons 701 and 703. Line 709, which appears thicker than line 707, indicates that the corrlation between neurons 703 and 705 is stronger than the correlation between neurons 701 and 703. The correlation map may be superimposed on an image array selected from the original input data. If the correlation coefficient is below a predetermined threshold amount, the corresponding line may be omitted from the correlation map. The predetermined threshold amount may be specified by a user or a default threshold may be used.

In still another exemplary embodiment, a cross correlogram is drawn to show potential causality among neuron activity. This can be used to find neurons with events that consistently precede or follow events of another neuron. A cross correlogram simply creates a histogram of the time intervals between events in two specified neurons. A line of height proportional to the number of times the first neuron is active one time interval following activity by the second neuron is plotted at +1 on the x-axis of the histogram. A line of height proportional to the number of times the first neuron is active two time intervals following activity by the second neuron is plotted at +2 on the x-axis of the histogram, and so on. An example of such a cross correlogram is illustrated in Fig. 6. The line 601 represents the number of occasions the first and second neurons were coactive, while line 607 represents the number of times the first neuron was active three time intervals after the second neuron was active. A cross correlogram may be performed on a single neuron to detect temporal characteristics in the neuron's firing such as the fact that the neuron is active with a period of every three time intervals a certain number of times during the observation.

IDL code implementing all of the preceding steps of the exemplary embodiment is attached hereto as Appendix A. The procedure "multicell" and "multicell_event" are the main procedures. All relevant sub-procedures and functions are also included in Appendix A.

In another exemplary embodiment of step 111 in Fig. 1, the data is analyzed by finding a hidden Markov state sequence from the second simplified data array. This embodiment uses the principal of Hidden Markov modeling described in Rabiner, A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proceedings of the IEEE, vol. 77 pp. 257-286 (1989). Essentially, a Markov model is away of modeling a series of observations as functions of a series of Markov states. Each Markov state has an associated probability function which determines the likelihood of moving from that state directly to any other state. Moreover there is an associated initial probability matrix which determines the likelihood the system will begin in any particular Markov state. In a hidden Markov Model, the Markov states are not directly observable. Instead, each state has an associated probability of producing a particular observable event. A complete Markov model requires the specification of the number of Markov states (N); the number of producible observations per state (M); the state transition probability matrix (A), where each element ay of A is the probability of moving directly from state i to state j; the observation probability distribution (B), where each element bj(k) of B is the probability of producing observation k while in state i; and the initial state distribution (P), where each element pj of P is the probability of beginning the Markov sequence in state i. In the exemplary embodiment, it is assumed that the number of times a neuron is active within each Markov state follows the Poisson distribution. Thus, each neuron in each state has an associated Poisson Lambda parameter, which can be understood in the exemplary embodiment to correspond to the neuron's average firing rate within the given Markov state. The set of all of these Lambda parameters is then assumed to be the B matrix. Given estimations of the Markov Model parameters, the method uses the Viterbi algorithm to find the single best state sequence, i.e. the state sequence that most likely occurred to generate the observed results. The parameter N may be selected by the user or a default number such as 4 states may be used. The Viterbi algorthim is described as follows:

Initialization

δ_λ(i) = p_ib_i(0_λ) l ≤ i ≤ N ,

where δ_t(i) is the highest probability along a single path at time t that accounts for the first t observations and ends in state i, and ψ is used to store the argument which maximizes δ_t(i). Recursion:

l ≤ j ≤ N,

Ψ, (/) = ^arS max[δ_t__r (t)α,_y ] 2 < t ≤ T l≤i≤N

l ≤ j ≤ N,

Termination:

αrg w ^_j. (Ϊ l≤i≤N

Path (backtracking):

Once a possible state sequence is generated, a cross-correlogram between neurons in a predetermined state can be plotted using the methodology previously described. The state may be selected by the user or a default state such as the first state may be used.

IDL code implementing the preceding embodiment involving the hidden Markov model is attached hereto as Appendix B. The procedure "hiddenjmarkov" and "hidden_markov_event" are the main procedures. All relevant sub-procedures and functions are also included in Appendix B.

In another exemplary embodiment the data is analyzed by performing a singular valued decomposition (SVD) on the data in the first simplified array formed in step 105. In this embodiment, steps 107 and 109 may be skipped. In a singular valued decomposition, the data set is reduced from N dimensions, where N is the number of selected neurons, to d dimensions, where d is specified and is less than N. The SVD algorithm, which is well known to one of ordinary skill in the art and is specified in the code in Appendix C, fits the observed data to a data model that is a linear combination of any d number of functions of the spaces of data (such as time and location). The SVD algorithm discards the eigenmodes corresponding to the smallest N-d eigenvalues. In this embodiment, the result that is plotted for visual analysis is the level of each neuron's contribution to each of the calculated d modes where a higher contribution to the mode is a larger number on the y-axis of the plot.

In a further exemplary embodiment not shown in Fig. 1, the SVD calculation described may be performed on the data in the first simplified data array before step 107. The data corresponding to the first eigenmode of the SVD, which will often correspond to noise in the original data, may then be removed from the first simplified data array. The modified first simplified data array, with noise removed, may then be processed beginning with step 107 in Fig. 1 as previously described. In this fashion, the data analysis, such as the HMM analysis could be performed on data with less noise, generating more useful results.

TDL code implementing the preceding embodiment involving the singular value decomposition algorithm is attached hereto as Appendix C. The procedure svd_gui and svd_gui_event are the main procedures. All relevant sub- procedures and functions are also included in Appendix C. Although the present invention has been described in detail, it should be understood that various changes, substitutions and alterations can be made hereto without departing from the scope or spirit of the invention as defined by the appended claims.

APPENDIX A

1997, (c) Copyright: Vik and Ted Sports Pvt. Ltd.

Name: cell_location.pro

Synopsis: cell_location

Description: Go into the directory of the movie you want to analyze .

; The program will find the location array and print it in a spreadsheet form with the cell numbers as well.

The first column is the cell numbers

The second column is the x value

The third column is the y value pro cell_location common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common choose_cell_convert_delta, filepathl cd, filepathl filel=pickfile (file= ' Cell_locations_exp_# ' /write) if (filel eq ' ' ) then begin mess=WIDGET_MΞSSAGE ( 'Saving locations cancelled. /INFORMATION) endif else begin no=n_elements (location. size) temp=intarr (3 , no) tem (0, *) =indgen(no) +1 temp (1, *) =location. coord (0) temp (2, *) =location. coord (1) close, 1 openw, 1, filel printf, 1, ¹ Cell Number printf, 1, temp close, 1 endelse end

Copyright (c) , 1997 by Vikram Sheel Kumar

NAME: choose_cells_chief

SYNOPSIS: choose_cells_chief

DESCRIPTION:

This is the program that creates the widget with which you can choose cells

,- It asks you to open the frame you want to pick cells from and then uses the filepath for that image as the root directory to find the remaining images

This is called by Multicell_event .

Modified on 7/13/97 by Vikram Sheel Kumar to add Find F/F button

Modified on 7/31/97 by VSK to put in check for first frame numbered with a 0 (taken out as only Movie.001 are left for many analysis) ; Modified by Rob, 2.99 to flip the tiff image that comes in, so it looks normal . ; Modified by Rob, 4.22.99, to generalize movie root name reading

pro choose_cells_chief common mother_com, pixel_array, yes_no_values, coef, location, cell_no , frame_no, fi1e_name , $ total_frame_no, time_resolution, x_size, y__size, box_size common test, str_ing, state3 common choose_cell_com, image, locations_defined common choose_cell_convert_delta, filepathl locations_defined = 0 name = dialog__j?ickfile (/read, title= ' For Choosing Cells, Select Movie Slice', get_j>ath=filepathl)

; ; determine movie root name- everything between filepath and the first 0 character is the root name. pos_first_zero = strpos(name, '0', strlen (filepathl) ) file_name = strmid(name, strlen (filepathl) , pos_first_zero - strlen (filepathl) )

WIDGET_CONTROL, state3.name, SET_VALUE=file_name if (name eq ' ' ) then begin mess = IDGET_MESSAGE ( ' No image was chosen', /INFORMATION) endif else begin image = tiff_read(name) ;this reads the first image to find the sizes image = reverse (image, 2) size_o_i = size (image) x_size = size_o_i(l) ,-this is the x size of the image y_size = size_o_i(2) ,-this is the y size of the image base3 = WIDGETJ3ASE (/COLUMN, title='Cell Choosing

Widget ' ) draws = WIDGET__DRAW (base3 , xsize=x_size, ysize=y_size) slide2 = WIDGET_SLIDER (base3 , value=3, maximum=60, minimum=l, title='Each box is 2x on a side') slide3 = WIDGET_SLIDER(base3, value=total_frame_no , maximum=total_frame_no, minimum=4, title= ' Total no of frames you shot ' ) buttonl = cw_bgroup (base3 , /row, ['Choose Cells', 'Find F/F ' ] , /return_name , UVALUE=3 ) WIDGET_CONTROL, /realize, base3 statel = {slide2 :slide2 , slide: slide3, draw:draw5} WIDGET_CONTROL, WIDGET_INFO (base3 , /Child), SΞT_UVALUE=Statel xmanager, ' choose_cells_chief ' , base3 endelse end

1997, (c) Copyright by Vikram Sheek Kumar NAME: choose_cells_chief_event SYNOPSIS: choose_cells_chief_event , event DESCRIPTION:

This is the event handler for the widget choose_cells_chief .pro ,- It allows the user to choose cells from a slice of a movie. It calls pick__box which is the program

; that tracks the mouse and finds where the cells are drawn. It also calles convert_to_delta_f_box which opens

; all of the frames from the movie and first finds the pixel values and then converts them into F/F values . MODIFICATION HISTORY

Created by Vikram Sheel Kumar in 1996 Modified: May 21, 1997 (graduation!) ; Added cell_location which saves cell locations to file you specify.

,- Modified: July, 13, 1997 : to allow you to change to file name for saving the cell locations .

; Modified: March 17, 1999: initialize threshold arrays for cell spiking, initialize yes__no_values , pixel_array, and coef arrays ; March 26, 1999: added ' original_data ' so we could go back' n¹ forth between filtered and unfiltered. ; May 10, 1999-. individual cell rms thresholds

EVENT HANDLER OF THE ABOVE pro choose_cells_chief_event, event common mother_com, pixel_array, yes_no_values, coef, location, cel1 no , frame_no, fi1e name , $ total_frame_no, time_resolution, x_size, y_size, box_size common old_skool_data, original_data common test, str_ing, state3 common with_choose_cells_com, thresholdl, threshold2, rms_threshold common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined common choose_cell_com, image, locations_defined common choose_cell_convert_delta, filepathl stateholder=WIDGΞT_INFO (event .handler, /Child) WIDGΞT_CONTRO , stateholder, GET_UVALUE=statel ,-find what is in /Child cell_no=0B frame no=0B CASE string (event. alue) OF 'Choose Cells' : BEGIN

WIDGET_CONTROL, statel.draw, GET_VALUE=win wset, win tvscl, image WIDGET_CONTROL, statel . slide2 , GET_VALUE=box_size

CD, filepathl location=pick_box (box_size, x_size, y_size) ; choose cells through this cell_no=location(0) .size ; location (0) since we only want one size cells_defined = 0 spikes_defined = 0 locations_defined = 1 mess=WIDGET_MESSAGE ( 'Press Choose Cells to reselect cells. Otherwise press Find F/F', /INFORMATION) END ' Find F/F ' : BEGIN if (locations_defined eg 0) then begin mess = WIDGET_MESSAGE ( 'You have to choose cells before I can find the F/Fs ! ' , /INFORMATION) endif else begin

WIDGET_CONTROL, statel . slide, GET_VALUE=frame_no cell_no = location (0) . size ; location (0) since we only want one size pixel_array = convert_to_delta_f_box (location, frame no, file_name)

WIDGET__CONTROL, event . to , /HOϋRGLASS cells_defined = 1 ,-this lets the programs know that cells have been defined! spikes_defined = 0 ,-this lets others know spikes not defined yet mess = WIDGET_MESSAGE (' Calculated the F/Fs. NOW saving the locations of the cells ... ' , /information) cell_location ,-this prints the cell locations to a file

; ; save original data for future restoration original_data = pixel_array

; ; create threshold arrays for each cell- we do it here just to initialize it as soon as cells are created

,- ; to get it out of the way rms_threshold = fltarr (cell no) thresholdl = fltarr (cell no) threshold2 = fltarr (cell no) for cell = 0, cell no - 1 do begin rms_threshold [cell] = 2.0 thresholdl [cell] = 2.0 threshold2 [cell] = 3.0 endfor ; ; init all other arrays that depend on values of frame_no and cell_no yes_no_values = intarr (cell_no, frame_no) coef = fltarr (cell_no, cell_no) ,- ; set new frame number if changed w/ slider WIDGET_CONTROL, state3. frame, SET_VALUE=frame_no endelse END

ELSE: BEGIN WIDGET_CONTROL, statel . slide2 , GΞT_VALUE=box_size

END ENDCASE end

NAME: choose_correl SYNOPSIS: choose__correl DESCRIPTION: ; This program creates the widget that calculates the cross correlation coefficients

; It is called by Multicell and calls choose_correl_event

; MODIFICATION HISTORY

; Created June, 1996 by VSK

WIDGET TO COMPUTE

CORRELATION pro choose_correl common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no , fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size common things_com, state3 base4=WIDGET_BASE (/ROW, title= ' Choose pvalue & windowsize¹) fieldl=CW_FIELD(base4, /RETURN_E EN S , /FLOATING, TITLE='p value for chi squared ' , $

VALUΞ=.05, UVALUE=0) field2=CW_FIELD(base4, /RETURN_EVΞNTS , /INTEGER, TITLE= ' window size, ONLY ODD¹ , $

VALUE=1, UVALUΞ=0) buttonl=WIDGET_BUTTON(base4, VALUE=' Go Cross Correlate', UVALUE=2)

WIDGET_CONTROL, /realize, base4 state3={fieldl:fieldl, field2 : field2, button :buttonl}

WIDGET CONTROL, WIDGET INFO (base4 , /CHILD), SET UVALUE=state3 xmanager, ' choose_correl ' , base4 end

; NAME:

; choose_correl_event, event ; SYNOPSIS:

; choose_correl_event, event

; DESCRIPTION:

; This handles the events for choose_correl . Together they allow the user to ,- find the cross correlation coefficients and specify a p value and bin size.

,- As a result, the cross correlation matrix is printed to a file in the working directory

7 ; MODIFICATION HISTORY

; Created - June, 1996 by VSK

; Last Modified - July 13, 1997 by VSK to allow the user to name the correlation coeff and

; write it to a specific directory and check to see if window is even or odd

EVENT HANDLER pro choose_correl_event, event ;__ global variables_ common mother_com, pixel__array, yes_no_values , coef , location, cell_no, frame_no, file_name, $ total_frame_no, time_resolution, x_size, y_size, box_size common things_com, state3 common flags , cells_def ined, spikes_defined, single_j>lot_defined, correl_coef_def ined stateholder3=WIDGET_INFO (event . handler, /CHILD) WIDGET_CONTROL, stateholder3 , GET_UVALUΞ=state3 WIDGET_CONTROL , state3 . f ieldl , GET_VALUE=set_ p_value

WIDGET CONTROL, state3 . f ield2 , GET VALUE=winsize

. calculate the correlation coefficients for all the cell pairs if ( (set_p_value It 0) OR (setj>_value gt 1) ) then begin mess=WIDGET_MESSAGE( 'The p must be between 0 and 1!', /ERROR) endif else begin if (winsize ge frame no) then begin mess=WIDGET_MESSAGE( 'The window size cannot be larger than the total number of frames in the movie ! ' , /ERROR) endif else begin if ( (winsize mod 2) eq 0) then begin mess=WIDGET_MESSAGE( 'Please enter an odd window size for uniformity¹, /INFORMATION) endif else begin coef=correl_coeff (yes_no_values, p=set_p_value, cells=cell_no, $ win=winsize, frame no) close, 1 openr, 1, 'correl.dat' stat_stuff=fstat (1) file_size=stat_stuff . size close, 1 print, file_size

;if (file_size gt 8112) then begin ;FIX FOR FILES

TOO LARGE LATER 12/5/97 Vsk

; mess=WIDGET_MESSAGE( 'File is too large to display through a widget. Open it manually. If you did not save it, it is named correl . dat ' , /INFORMATION)

; correl_coef_defined=l ; endif else begin xdisplayfile, 'correl.dat', title = "Correlation Coefficient Matrix", group = event. top, width = 85, height = 45 correl_coef_defined=l file2=pickfile (/write, file= ' Correlation_Coeff_exp_#_' ) if (file2 eq ' ' ) then begin mess=WIDGET_MESSAGE ( 'The Correlation Coefficient Matrix has not been saved. ' , /INFORMATION) endif else begin get_lun, lunl openw, lunl, file2 printf, lunl, coef free_lun, lunl close, lunl

;endelse endelse endelse endelse endelse end pro choose_crosscorr base=Widget_base (/column) buttonl=widget_button(base, value= ' Standard Cross Correlogram¹ ,uvalue=l) button2=widget_button(base, value='Cell firing rate Cross Correlogram' ,uvalue=2) widget_control,base, /realize xmanager, ' choose_crosscorr ' , base end pro choose_crosscorr_event, iv widget_control , ev.id, get_uvalue=uval case uval of

1 .- draw_cross 2 : draw_cross2 endcase end

NAME: choose_threshold

SYNOPSIS : choose threshold

DESCRIPTION:

This program creates the widget with which the user can change the thresholds for the

; detection of spikes in cells. It is called by calling Find Spikes to Multicell.

; It uses MAKE_BINARY (see choose_threshold_event) to make the binary array .

MODIFICATION HISTORY:

Created on October 3, 1996 by Vikram Sheel Kumar Modified on Jul 13, 1997 by Vikram Sheel Kumar to complete documentation .

Modified, 2.99, by Rob to make threshold vars into arrays

WIDGET TO COMPUTE

CORRELATION pro choose_threshold

THE GLOBAL VARIABLES common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame__no, file name, $ total_frame no, time_resolution, x_size, y_size, box_size common with_choose_cells_com, thresholdl, threshold2, rms_threshold common things_com, state3 base5=WIDGΞT_BASE (/Column, title= ' Choose threshold for spike detection' ) fieldl=CW_FIΞLD(base5, /RETURN_ΞVENTS, /FLOATING, TITLE= ' Threshold for 2 consec. frames to be a spike ',$ VALUE=threshl, UVALUE=0) field2=CW_FIELD(base5, /RΞTURN_EVENTS , /FLOATING, TITLE= ' Threshold for 3 consec. frames to be a spike ', $ VALUE=thresh2 , UVALUE=0) buttonl=WIDGET_BUTTON (base5 , VALUE= ' Find Spikes ' , UVALUE=2 )

WIDGET_CONTROL, /realize, bases state3={fieldl:fieldl, field2 :field2, butto :buttonl}

WIDGET_CONTROL, WIDGET_INFO (base5 , /CHILD), SΞT_UVALUE=state3 xmanager, 'choose_threshold' , bases

END

NAME: choose_thresho1d_event

SYNOPSIS: choose_threshold_event (called by xmanager)

DESCRIPTION:

This is the event handler for the choose_threshold widget.

With this, the user can change the threshold value used in determining which increases in

Calcium intensity correspond to action potentials.lt uses make_binary .pro to convert the pixel or deltaF arrays into 0 or Is .

MODIFICATION HISTORY

Created on October 3. 1996 by Vikram Sheel Kumar

Modified on Oct. 31, 1996 by VSK to incorporate deltaF values

Modified on 2.16.97 by VSK to save data to a file named by the experiment number

Modified on 3.12.97 by VSK to be able to calculate different thresholds for different cells.

Modified on 7.13.97 by VSK to set tests to ensure stability of working

Modified, 2.99, by Rob to make threshold vars into arrays

EVENT HANDLER pro choose_threshold_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common things_com, state3 common with_choose_cells_com, thresholdl, threshold2, rms_threshold common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined stateholder3=WIDGET_INF0 (event .handler, /CHILD) WIDGET_C0NTR0L, stateholder3 , GET UVALUE=state3 WIDGET_CONTRO , state3. fieldl, GET_VALUE=threshl WIDGET_CONTROL, state3. field2 , GΞT_VALUE=thresh2 if ( (threshl eq 0) OR (thresh.2 eq 0)) then begin mess=WIDGET_MΞSSAGE ( 'Threshold cannot be 0!¹, /Error) endif else begin

WIDGET_CONTROL, event . top, /hourglass for cell = 0, cell_no - 1 do begin thresholdl [cell] = threshl threshold2 [cell] = thresh2 endfor yes no_values=make_binar (pixel_array, frame no, cell_no, thresholdl, threshold2) spikes_defined=l ,-this lets other programs now that spikes have been defined mess=WIDGET_MESSAGE( 'Done finding spikes boss.', /INFORMATION) endelse end ; COPYRIGHT (C) 1997-1997 BY VIKRAM SHEEL KUMAR

NAME: contingency_values

SYNOPSIS: contingency_values, filearrayl, filearray2

DESCRIPTION:

This program creates a contingency table for the data of two spike trains and from it gives a chi squared value.

METHOD: observed_array is an array of 4 elements they are as follows :

; observed_array (0) =number of hit_hits, i.e the number of times the first cell and the second cell have spikes at exactly the same time. observed_array (1) = number of observed hit_misses observed_array (2) = number of observed miss_hits observed_array (3) = number of observed miss_misses

This is called by CORREL_COEFF . PRO to calculate the correlation coefficients between all the cells of

,- an analysis. Using the contingency table, we get a chi squared value which we compare with the chi squared

; value created when the user specifies the p for the data to be significant. Each pair of cells has its ; own contingency value . /

; MODIFICATION HISTORY:

; Created by Vikram Sheel Kumar in May, 1996

; LAST MODIFIED by Vikram Sheel Kumar on 7-13-97 to complete documentation

function contingency_values, filearrayl, filearray2 ,-this program will print out a contingency table and then print out the chi-square value.

;it also print out the probability that the null hypothesis is true ,-Null hypothesis: The two files are independant

/Constraint: THIS WILL ONLY WORK IF THE SPIKE TRAINS FOR THE CELLS ARE IN TERMS OF 0 AND 1

_/declaring the arrays observed_array=fltarr (4) expected_array=fltarr (4) result=fltarr (4)

_/This makes a spike for celll = 2 where a spike for cell2=l. Therefore celll-cell2 will be 0 only ,-at times when both celll and cell2 have a 0 (at miss_miss) . A

_/Due to this, a hit_miss can be found by seing where the above difference array is equal to 1.

_/eg., 1 0 0 1 0 1 -> 2 0 0 2 0 2 filearrayl=2*filearrayl

_/the following array is used to calculate all of the values needed in the table

_/we subtract the spike trains of the two cells into difference_array.

_/therefore if difference_array has a 0, there difference_array=filearrayl-filearray2

_/this is the m(l,l) values or the hit_hit

,-a spike in filearrayl will be a 2. That in filearray2 will be a 1. Therefore, in the difference array ,-a hit_h.it will be represented by a (2-1) or a 1 arr_hit_hit=difference_array eq 1 if (total (arr_hit_hit) gt 0) then begin observed_array (0) =total (arr_hit_h.it) endif else begin observed_array (0) =0 endelse

,-this is the m(0,0) values or miss_miss arr_temp3=difference_array eq 0 if (total (arr_temp3) eq 0) then begin observed_array (3) =0 endif else begin miss_miss_data=total (arr_temp3) observed_array (3) =miss_miss_data endelse

,-observed__array (3) =n_elements (difference_array (where (difference _array eq 0) ) )

/this is the m(l,0) values or the hit_miss, i.e the number of times when at an instant in time, /the first cell has a spike, but the second cell does not. arr_hit_miss=difference_array eq 2 if (total (arr_hit_miss) eq 0 ) then begin observed_array (1) =0 endif else begin hit_miss_data=total (arr_hit_miss) observed_array (1) =hit_miss_data endelse ,-this is the m(0,l) values or the miss_b.it, i.e the number of times when at an instant in time, the

/second cell has a spike while the first cell does not. (celll (t)=0, but cell2(t)=l) arr_miss__hit=difference_array It 0 if (total (arr_miss_hit) eq 0) then begin observed_array (2) =0 endif else begin miss_hit_data=total (arr_miss_hit) observed_array (2) =miss_hit_data endelse sum_first_row=observed_array (0) +observed_array (1) sum_second_row=observed_array (2) +observed_array (3) sum_first_col=observed_array (0) +observed_array (2) sum_second_col=observed_array (1) +observed_array (3 ) n_total=sum_first_row+sum_second_row

_/now we calculate the observed and expected values for chi squared

/for j=l, i=l expected_array (0) = (sum_first_row*sum_first_col) /n_total ,-for i=l, j=0 expected_array (1) = (sum_first_row*sum_second_col) /n_total _/for i=0, j=l expected_array (2) = (sum_second_row*sum_first_col) /n_total _/for i=0, j=0 expected_array(3) = (sum_second_row*sum_second_col) /n__total for k=0, 3 do begin result (k) = ( (observed_array (k) - expected_array (k) ) * (observed_array (k) - expected_array (k) ) ) /expected_array (k) endfor answer=total (result) return, answer end

NAME: convert_to_delta_f_box

SYNOPSIS: convert_to_delta_f_box, cell_location, frame_no, filename

DESCRIPTION: _/ This program converts the pixel values into deltaF/F values to normalize the data. The deltaF/F value varies from 0 to 100%. By using this technique of measuring intensity, _/ a brighter cell will need a greater pixel change than a dimmer cell to achieve the same

_/ deltaF/F. The baseline is created by averaging the first 4 frames . _/ Formula used:

Delta F F - F(0) _{= 10}o

F F ( 0) ;

_/ MODIFICATION HISTORY:

; Created by Vikram Sheel Kumar in November, 1996 Modified on 12/3/96

_/ Modified on 7/13/97 to complete documentation ,- Modified by Rob, 2.99 to flip the tiff image so that it's read correctly

; Modified by Rob, 3.16.99: deleted the delta_f_main common block: had one unshared variable, 'baseline'...

,- as far as I can tell, it isn't used anywhere else, so: see ya ! function convert_to_delta_f_box, cell_location, frame_no, filename2 common choose_cell_convert_delta, filepathl

!except=2 realname=strarr (4 ) realname (0) =filename2+string( ' 00 ' ) _/if the frame_number is between 1 and 9 the units place of it's log will 0 and this will be called realname (1) =filename2+string( ' 0 ' ) _/if the frame_number is between 10 and 99, byte (it's log) =1 and this will be called realname (2) =filename2 _/if the frame_number is between

100 and 999, byte (it's log) =2 and this will be called realname (3 ) =filename2 cell_no=cell_location(0) .size

/incorporate all this into a structure pixel_values=fltarr (cell_no, frame_no) delta_f_values=fltarr (cell_no, frame_no) half_side=cell_location(0) .half_side _/HALF LENGTH OF BOX SIDE total_side=half_side*2+l _/FOR EG.

HALF_SIZE=3, TOTAL PIXELS CALCULATED OVER =7*7=49 area=total side^{^}2

OPEN EACH FRAME AND MEASURE THE PIXEL

VALUE OF EACH CELL THERE for frame_number=0, frame_no-l do begin tempo=tiff_read(strcompress (realname (byte (aloglO (frame_number) ) ) +string (frame_number) , /remove_all) ) tempo=reverse (tempo, 2) for cell_numb=0, cell_no-l do begin x=cell_location (cell_numb) . coord ( 0 ) y=cell_location (cell_numb) . coord ( 1) box_array=tempo (x-half_side :x+half_side , y- half_side : y+half_side) pixel_values (cell_numb , f rame_number) = (total (box_array) ) /area _/ this is the average pixel value in a box endfor endfor

CONVERT EACH PIXEL ARRAY INTO AN ARRAY OF

DELTA F VALUES for cell_number = 0 , cell_no - 1 do begin baseline = total (pixel_values (cell_number, 0:3)) / 4.00 _/average first 4 frames to get start delta_f_values (cell_number, *) = (pixel_values (cell_number, *) - baseline) / (baseline * 1.0) endfor delta_f_values (*, *) =TEMPORARY(delta_f_values (*, *))*100 return, delta_f_values end

NAME: correl_coeff .pro

SYNOPSIS: correl_coeff, yes_no_values, P_VALUΞ=p_value, CELL NUMBER=cell number, WINSIZE=winsize

_/ INTRODUCTIO :

_/ This routine is used in the widget MULTICELL to output a correlation coeff . matrix with

_/ the insignificant correlation coeff =0. It uses the functions Contingency_values and find_matches from

,- the directory MULTICELL. This procedure is called by the widget CHOOSE__CORREL .

7

_/ INPUTS: ,- p_value-> This is the probability that a random variable is greater than the cut off value in

,- a chi squared distribution. For example, p=.05 means the cut off value should be high enough

_/ so there is only a 5 percent chance that any random number can be greater.

_/ cell_number-> the number of cells chosen from the slice.

_/ winsize-> the size of the window to look for coincident spike in. The default is 1

DESCRIPTION:

The function does the following:

1. Uses the IDL function chi_sqrcvf to find the cut off value for the input p value. It is _/ called as chi_sqrcvf (p_value, degrees of freedom) / 2. Creates a contingency table for a pair of cells by the home-made function

,- contingency_values . pro. From this table, it then finds the chi squared value for the null hypothesis that

_/ the cells are independant.

,- 3. It prints out all of the correlation coefficients that are significant (fail the null

/ hypothesis, but pass the cut_off test) and indicate the cells of the pair are dependant .

OUPUT : a matrix of the filtered correlation coefficients .

MODIFICATION HISTORY:

Created in July, 1996 by Vikram Sheel Kumar Modified on 7/8/96 by Vikram Sheel Kumar Modified on 7/13/97 by Vikram Sheel Kumar to update documentation ,- Last Modified on 7/18/97 by VSK to allow the contingency table to be made with a binsize other than 1 . MAHj PROGRAM

_/NOTE: To see all of the cross correlation coefficients, make the p value=l! function correl_coeff , yes_no_values, p=p_value, cells=cell_number, win=winsize, frame no forward_function contingency_values, find_matches output__array=fltarr (cell_number, cell_number)

/ FIND THE CUT OFF VALUE-

limit=chisqr_cvf (p_value, 1) ,-the degree of freedom is 1 in our case

. CONTINGENCY TEST AND CORREL COEFF

CALC . bin_no=frame_no/winsize _/using the window size, this is the number of bins left_overs=frame_no mod winsize _/these are the frames that didn't fit in the bins, but are still in the array (at the end) total_elements=bin_no+left_overs _/this is the total number of elements in the binned array if (total_elements It 1) then begin mess=WIDGET_MESSAGE ( 'Window size too large or two few frames.',

/ERROR) endif else begin close, 1 IF (winsize eq 1) THEN BEGIN _/IF THE WINDOW

SIZE IS ONE, NO BINNING IS DONE! for i=0, cell_number-l do begin for j=0, cell_number-l do begin temp__value=contingency_values (yes_no_values (i, *) , yes_no_values (j , *) ) if (temp_value ge limit) then begin output_array (i , j ) =find_matcb.es (yes_no_values (i, *) ,yes_no_values (j , *),$ winsize) endif endfor endfor ENDIF ELSE BEGIN

celll_binned=intarr (total_elements) cell2_binned=intarr (total_elements) stop_index=-l index=0

,-in case no bins are made, then the stop index used in the second for loop below will be used for i=0, cell_number-l do begin for j=0, cell_number-l do begin for bin_number=0, bin_no-l do begin _/making the binned cell array start__index=bin_number*winsize _/this is the starting point for a bin stop_index=start_index+winsize-l

_/this is the ending point for a bin celll_binned(bin_number) =total (yes_no_values (i, start_index:stop _index) ) ,-we are binning the input arrays for the contingency table cell2_binned (bin_number) =total (yes__no_values (j , start_index: stop _index) ) ,- for the second cell being looked at endfor index=0 for bin_number=stop_index+l, frame_no-l do begin

,-add the elements that don't fit in the bins celll_binned (bin_no+index) =yes_no_values (i , bin_number) cell2_binned(bin_no+index) =yes_no_values (j , bin_number) index=index÷1 endfor celll_binned=celll_binned ge 1 ,-this converts the arrays into binary ones (eg if three elements in a bin each have a value 1 then the total will be 3 that will be converted to 1 here) cell2_binned=cell2_binned ge 1 temp_value=contingency_values (celll_binned, cell2_binned) if (temp_value ge limit) then begin output_array (i , j ) =find_matches (yes_no_values (i, *) ,yes_no_values (j , *) , $ wmsize) endif endfor endfor

ENDELSE close, 1 openw, 1, 'correl.dat' / CHOOSE_CORREL USES XDISPLAYFILE AND THIS FILE TO WRITE printf, 1, output_array close, 1 return, output_array ,- this is the correl coeff matrix endelse end

Unauthorized reproduction prohibited.

NAME:

CORREL MAP IMAGE

PURPOSE :

This procedure creates a correlation map between the input cells. A slice of the movie you are analyzing is loaded and the cells you chose that are correlated in either direction are joined by lines whose thicknesses are proportional to their correlation coefficients

MODIFICATION HISTORY:

Created in July, 1996 by VSK Modified on Sept 25, 1996 by VSK

Modified on Oct 10, 1996 by VSK to changed tv into tvscl Modified on Jan 29, 1997 by VSK so names of cells can be read. Modified on July 11, 1997 to allow user to choose frame on which to plot map and make sure the chosen slice is from the correct movie.

Modified on July, 13, 1997 by VSK to update documentation

Modified 2.99, Rob, to flip the tiffs!

Modified, 4.05.99, Rob, to move the numbering pro correl_map_image common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size name_of_file=strcompress (file_name (0) -i-string ('*')) name=pickfile (FILE=file_name (0) , FILTER=name_of_file) print, name if (name eq ' ') then begin _/if you did not select an image mess=WIDGΞT_MΞSSAGE{ 'No image was selected¹, /INFORMATION) endif else begin image=tiff_read (name) image=reverse (image, 2)

_{/ /} modest error checking- this doesn't actually check to see that the movie is the

; ; same as first picked, only that the sizes agree. size_info=size (image) x=size_info (1) y=size_info (2) if ( (x ne x_size) OR (y ne y_size) ) then begin mess=WIDGET_MESSAGE( 'Wrong movie! With Multicell, however, you do get a second chance.') endif else begin loadct, 0 base4=WIDGET_BASE (/COLUMN, title= ' Correlation Map on Image ' ) draw6=WIDGET_DRAW(base4, xsize=x size, ysize=y_size)

WIDGET_CONTROL, /realize, base4 tvscl, image for s=0, cell_no-l do begin x=location(s) .coord (0) _/the x axis of where you clicked on cell y=location(s) .coord (1) _/the y axis of where you clicked on cell xyouts, x-50, y, s+1, alignment=.5, /device for t=0, cell_no-l do begin if (coef (s, t) gt 0) then begin plots, [location (s) . coord (0) , location (t) . coord (0) ] , [location (s) . coord (1) , location (t) .coord (1) ] , $ thick=3*coef (s, t) , /device endif endfor endfor endelse endelse

!p.font=0 end

7

,- NAME:

CORREL_MAP_IMAGE_PLANE _/ PURPOSE:

_/ This procedure creates a correlation map between the input cells. A correlation i,

_/ a circular representation of correlation between cells of the slice. Those cells that _/ are correlated in either direction are joined by lines which are proportional to their

_/ correlation coefficients PARAMETERS : number: Number of cells that have been chosen from the cell. MODIFICATION HISTORY:

Modified: 2/16/97 by VSK to put the numbering back in. Modified: 7/10/97 by VSK to only draw dots for the active cells .

_/ Modified: 7/13/97 by VSK to update documentation _/ Modified: 4/05/99, Rob, to move the numbering pro correl_map_image_plane, coef_array common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size _//common with_correl_map_plan, cluster_array for s=0, cell no-1 do begin

_/ ,- 1) draw boxes over active cells

,- ,- we only worry about the active cells- cells that have at least one spike if total (yes_no_values (s, *) ) gt 0 then begin x=location(s) .coord (0) _/the x axis of where you clicked on cell y=location(s) . coord (1) _/the y axis of where you clicked on cell for j=-box_size, box_size do begin plots, [x+j , x+j] , [y-box_size, y+box_size] ,

/device endfor xyouts, x-50, y, s+1, font=-l, charsize=l.9, /device

_{/ /} 2) draw lines between 'correlated' cells for t=0, cell_no-l do begin if (coef_array (s, t) gt 0) then begin plots, [location (s) . coord (0) , location (t) . coord (0) ] , [location (s) . coord (1) , location (t) .coord (1) ] , $ thick=2*coef_array (s, t) , /device endif endfor endif endfor !p.font=0 end

(c) 1997, Copyright by Vikram Sheel Kumar Name: count_random_hits_2_manyX, least_no_of_matcb.es Synopsis : count_random__hits_2_manyX, random_array, least_no_of_matches _/ Description: This program tests to see how many times more than least_no_of_matches two random cells fire simultaneously. Description of variables : yes_no_values : the binary data from the actual slice _/ random_array : a random array created by RANDOM_TEST looking at the spikes in yes_no_values ,- least_no__of_matches : Minimum number of times that two cells must fire to be considered a coactive pair

; window_size: number of frames before and after to look for a coincident spike in the other cell ,- Created by VSK on June 1, 1997

,- Modified by VSK on June 3, 1997 so that it only counts the number of times pairs of cells of the

_/ population fire together more than or equal to least_no_of_matches times . ; Modified, 3.27.99, Rob, to add window_size and supporting code... ,- Modified, 4.12.99, Rob, to clean up code and improve speed ,- Modified, 4.24.99, Rob, speedup to ignore silent cells _/ Modified, 5.31.99, Rob, shifting to avoid wrapping function count_random_hits__2_manyX, random_array, least_no_of_matcb.es, window_size, num_active_cells common mother_com, pixel_array, yes_no__values, coef, location, cell_no , frame no, fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size random_hits=0 for cell_l_counter = 0, num_active_cells - 2 do begin for cell_2_counter = cell_l_counter + 1, num_active_cells

- 1 do begin temp_hits = 0 for window_counter = -window_size, window_size do begin temp_cell_l = intarr (frame_no + (2 * window_size) ) temp_cell_2 = intarr (frame_no + (2 * window size) ) temp_cell_l (window_size : frame_no+window_size- 1) = random_array (cell_l_counter, *) temp_cell_2 (window_size : frame_no+window_size- 1) = random_array (cell_2_counter, *) temp_cell = (temp_cell_l * shift (temp_cell_2, window_counter) ) temp_hits = temp_hits + total (temp_cell) endfor if (temp_hits ge least_no_of_matches) then begin random_h.its = random_hits + 1 endif endfor endfor return, random_hits end

1997, (c) Copyright by VSKaze et al. Inc.

Name : count_random_matches

Synopsis: count_random_matcb.es, random_array, times_repeated

Description: this program takes in a random array of spike trains and counts the number of sets of coactive spikes

; (where the number of spikes in a set is ' times_repeated' , e.g. if 2, we count the number of all possible pairs of coactive spikes

Defintion of variables : random_array is the array of the location of spikes for the simulated data

_/ times_repeated is the number of times the spikes are clustered together in calculating the

_/ combination. For example, if looking at pairs of two, times_repeated=2. If looking for cells firing 3X together, it is 3.

Formula Used: C(n,r)= n!/(r! * (n-r) ! )

Created by Vikram Sheel Kumar on May 22, 1997 ,- Modified, Rob, 1.00, windowing for gen sig test function count_random_matches , random_array, times_repeated, window_size common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size

,-find the number of hits for the simulated data. _/ we can do this by summing the rows (cells) of one column (frame) since spikes are

,- represented by 1 and then finding the combinations of this taken times_repeated at a time

_/ eg. cell 1 : 001100110101010100101010 _/ cell 2: 001000100001000001000001

_/ If in frame 23, only cells 2, 3, and 8 fire, the sum of the frames will

_/ be 3 and then the combination taken two at a time is C(3,2)=3 random_match.es = 0 if (window_size eq 0) then begin for frames = 0, frame_no - 1 do begin matches = total (random_array (*, frames)) _/finds n or the number of hits at this frame if ( (matches gt 0) and (matches ge times_repeated) ) then begin

_/ ,- find number of possible pairs matches_pair = factorial (matches) / (factorial (times_repeated) * factorial (matches-times_repeated) ) endif else begin matches_pair = 0 _/cant find factorial of a negative number endelse random_matcb.es = random_matches + matches_pair endfor endif else begin

_/ / first compress spike trains according to window size binned_random_array = intarr (n_elements (random_arra [*, o] ) , frame_no - (window_size * 2)) for frame_counter = window_size, frame_no - (window_size + 1) do begin for window_counter = -window_size, window_size do begin for cell_counter = 0 , n_elements (random_array [*, 0] ) - 1 do begin if {random_array [cell_counter, frame_counter + window_counter] ge 1) then begin binned_random_array [cell_counter, frame_counter - window_size] = 1 endif endfor endfor endfor

,- ; then look for matches in the binned spike trains (although a single spike can fall into many bins) for frame_counter = 0, n_elements (binned_random_array [0,

*] ) - 1 do begin matches = total (binned_random_array (*, frame_counter) ) _/finds number of hits at this frame if (matches ge times_repeated) then begin matches_pair = factorial (matches) /

(factorial.(times_repeated) * factorial (matches - times_repeated) ) ; finds number of possible pairs endif else begin matches_pair = 0 ,-cant find factorial of a negative number endelse random_matches = random_matches + matches_pair _/increases the matches count for the random data by that of this frame endfor endelse return, random_matches end

1997, (c) Copyright by VSKaze et al. Inc.

Name : count_random_matches_many

Synopsis : count_random_matches_many, random_array, times_repeated

Description: this program takes in the actual binary spike train and the random one and sees how ,- the actual one compares to the random one. To test the significance, it looks at the number of

_/ coincident spikes boths sets of data have and creates a uniform distribution for the random data _/ to see where the actual data lies. ,- Defintion of variables: yes_no_values is the binary spike train of the true data

/ random_array is the array of the location of spikes for the simulated data

Created by Vikram Sheel Kumar on May 22, 1997

Modified by Rob, 4.20.99 to change comparison from ' ge ' to what the user wants- either 'ge' or 'eq'

Modified, Rob, 5.8.99, window code

Modified, Rob, 5.14.99, binning for window code

Modified, Rob, 5.26.99, revision of binning function count_random_matches_many, random_array, times_repeated, ge_or_eq_test, window_size common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size random_matches = 0 if (window_size eq 0) then begin for frames = 0, frame_no - 1 do begin matches = total (random_array (* , frames)) _/finds number of hits at this frame if ( ( (ge_or_eq_test eq 0) and (matches ge times_repeated) ) or ( (ge_or_eq_test eq 1) and (matches eq times_repeated) ) ) then begin random__matches = random_matches + 1 endif endfor endif else begin

_/ / first compress spike trains according to window size binned_random_array = intarr (n_elements (random_array [* , 0] ) , frame_no - (window_size * 2)) for frame_counter = window_size, frame_no - (window_size + 1) do begin for window_counter = -window_size, window_size do begin for cell_counter = 0, n_elements (random_array [*, 0] ) - 1 do begin if (random_array [cell_counter, fratne_counter + window_counter] ge 1) then begin binned_random_array [cell_counter, frame_counter - window_size] = 1 endif endfor endfor endfor

,- ,- then look for matches in the binned spike trains (although a single spike can fall into many bins) for frame_counter = 0, n_elements (binned_random_array [0,

*] ) - 1 do begin matches = total (binned_random_array (*, frame_counter) ) /finds number of hits at this frame if ( ( (ge_or_eq_test eq 0) and (matches ge times_repeated) ) or ( (ge_or_eq_test eq 1) and (matches eq times_repeated) ) ) then begin random_match.es = random_match.es + 1 endif endfor endelse return, random_matches end

Name : delete_spikes_widget

Description: allows me to delete from all cells a spike in a particular frame . Basically just a time-saving device.

History: Created 3.1.99 by Rob.

Global units threshold for spike deletion added, 3.26.99, Rob. Delete all trailing spikes control added 4.1.99, Rob. pro delete_spikes_widget common delete_spikes, statel2 basel2=WIDGET_BASE (/Column, title= ' Spike Deletion') buttonll=WIDGΞT_BUTTO (basel2 , VALUE= 'Delete All Trailing Spikes', UVALUE=11) fieldl2=CW_FIELD(basel2, /RETURN_E ENTS , /INTEGER, TITLE= ' Frame for Global Spike Deletion', VALUE=0, UVALUE=0) buttonl2=WIDGΞT_BUTTON(basel2, VALUE= ' Delete Spike in This Frame ' , UVALUE=2 ) fieldl3=CW_FIELD(basel2, /RETURN_EVΞNTS , /FLOATING, TITLΞ= ' Threshold in dF/F Units for Global Spike Deletion', VALUΞ=0.0, UVALUE=0) buttonl3=WIDGET_BUTTON(basel2, VALUE= 'Delete Spikes Less Than Threshold', UVALUE=2)

WIDGΞT_CONTROL, /realize, basel2 statel2={buttonll:buttonll, fieldl2 :fieldl2, buttonl2 :buttonl2, fieldl3 : fieldl3 , buttonl3 :buttonl3 } WIDGET_CONTROL, WIDGET_INFO (basel2 , /CHILD), SET_UVALUΞ=Statel2 xmanager, ' delete_spikes_widget ' , basel2 end

Name : delete_spikes_widget_event

Description: allows me to delete from all cells a spike in a particular frame .

; Basically just a time-saving device. Obviously this is the event handler for the

' delete_spikes ' widget . Crude but it works .

History: Created 3.1.99 by Rob.

Global units threshold for spike deletion added, 3.26.99, Rob. Delete all trailing spikes control added, 4.1.99, Rob.

pro delete_spikes_widget_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common delete_spikes, statel2 stateholderl2=WIDGET_INFO (event .handler, /CHILD) WIDGET_CONTROL, stateholderl2 , GET_UVALUE=statel2 WIDGET_CONTROL, statel2. fieldl2 , GET_VALUE=frame_to_kill

WIDGET_CONTROL, statel2. fieldl3 , GET_VALUE=spike_unit_threshold

_{/ /} deleting ALL trailing spikes in ALL cells if (event. id eq statel2.buttonll) then begin for cell_counter = 0, (cell_no - 1) do begin for frame_counter = (frame_no - 1) , 2, -1 do begin if ( (yes_no_values [cell_counter, frame_counter] eq 1) and ( (yes_no_values [cell_counter, frame__counter - 1] eq 1) or (yes_no_values [cell_counter, frame_counter - 2] eq 1) ) ) then begin yes_no_values [cell_counter, frame_counter] = 0 endif endfor endfor endif

,- ,- delete spike in every cell that is in a given frame if (event. id eq statel2.buttonl2) then begin if (frame_to_kill ge 0) and (frame_to_kill It frame_no) then begin for cell_counter = 0, (cell_no - 1) do begin yes_no_values [cell_counter, frame_to_kill] = 0 endfor endif endif

,^• ,- delete all spikes less than a raw dF/F threshold of intensity if (event. id eq statel2.buttonl3) then begin for cell_counter = 0, (cell_no - 1) do begin for frame_counter = 0, (frame_no - 1) do begin if (pixel_array [cell_counter, frame_counter] It spike_unit_threshold) then begin yes_no_values [cell_counter, frame_counter] = 0 endif endfor endfor endif end

Name : draw_3D_plot

Description: draws all frames of all cells, with intensity of dF/F on the z-axis,- i.e. it makes a pretty cool plot !

History: Created February 24, 1999, by Rob. pro draw_3D_plot, pixel_array window, 4, title='All Cells-3D' surface, transpose (pixel_array) , /horizontal end

Copyright (c) Synapse Innovation Inc. NAME : draw_cgram2

SYNOPSIS: draw_cgram2 , pixel_array, file_numberl, file_number2 , bin_size, cell_no

_/ INTRODUCTION: / This procedure is used to draw cross correlograms between any two cells. A cross correlogram

,- is a histogram of the time intervals between the two cells . If there is a peak at 0, the cells are ,- very highly correlated. This draws a cross correlogram between cell 2 with respect to cell 1. That is ,- a peak at +1 means that cell 2 leads cell 1 by 1 time frame.

/ INPUTS : ,- pixel_array. This is the array created in either pixel_binary_cells or pixel_binary_box

,- depending on the choosing method you used. It contains the pixel values of the cells and their

,- converted binary form. It is the concated array of: ;

7

; pix_array=intarr (-cell no, frame no) {actually called pixel_array in those pro} yes_no_values=intarr (cell_no, frame_no)

_/ where cell_no is the number of cells you chose

; and frame_no is the number of frames you chose to analyze

_/ therefore, what we are talking about is:

7 ; pixel_array =intarr (cell_no+cell+_no, frame_no)

7

; to find the data on cell x,

; pixel_array (x-1, *) will give the pixel values over time (from pix_values) ,- pixel_array (cell_no+x-l, *) gives 1 for a peak else

0 (from yes_no_values)

7

,- cell_numberl : This is the number of the first cell you want to cross correlate from ,- cell_number2 : This is the number of the second cell you are considering

_/ binsize: This is the window size you are calculating with. If any spikes fall in this, they are

_/ considered coincident. Default is 1 _/ cell_no: Total number of cells you have chosen

_/ ime: The time between two frames shot by the CCD camera (input from widget)

7

_/ DESCRIPTION: / This program looks at the yes_no_values part of pixel_array

(look above) and converts the

,- binary array into an array holding the positions of the spikes . It uses the in built function

; where () to do this. Then you use the function FIND_DISTANCE to calculate the time differences between

/ the spikes. It's source code is located in the same directory as this, 'All you need'. Briefly what

/ it does is it subtracts the two time arrays

7 / POSSIBLE IMPROVEMENTS: This doesn't have a constant max y-axis value for the graph. So if the max / is 1, all you see are straight verticals line. MODIFICATION HISTORY:

December, 19 1996 by Vikram Sheel Kumar Modified:

February, 27 1997 by VSK to allow one to change the scale of the xaxis when plotting _/ Last Modified :

; July, 10, 1997 by VSK to fix the bug that occured when the cells had no spikes

AUTHOR: Vikram Sheel Kumar vsk6@columbia.edu vikram@cubsps .bio . Columbia . edu

pro DRAW_CGRAM, yes_no_values , cell_numberl, cell_number2, bin_size, cell_no, time, xmin, xmax

,-this procedure will make a correlogram between the two files input. ,-the functions it uses are :findgen (size) *bin_size+mintemp ,-openfile, find_distance

!p.multi=0 forward_function FIND_DISTANCE write=string(format= ' ("cross-correlogram of cell",/, (iO) ) ', cell_numberl) write2=string(format= ' (" vs cell ",/,(i0))', cell_number2) timel=where (yes_no_values (cell_numberl-l, *) gt 0) ,-this finds where the spikes are time2=where (yes_no_values (cell_number2-l, *) gt 0) IF ((timel(O) eq -1) OR (time2 (0) eq -1)) then begin mess=WIDGET_MESSAGΞ ( 'One of the cells has no spikes') ENDIF ELSE BEGIN distemp=FIND_DISTANCE (timel, time2) mintemp=min (distemp) maxtemp=max (distemp) xmin=float (xmin) /time xmax=float (xmax) /time size=xmax-xmin IF bin_size It 1 THEN BEGIN mess=WIDGET_MESSAGE( 'Binsize must be >=1')

ENDIF ELSE BEGIN hist=histogram (distemp, binsize=bin_size, min=xmin, max=xmax) plot, (findgen(size) *bin_size+xmin) *time, hist, psym=10, title=write (0) ÷write (1) +write2 (0) +write2 (1) , $ xtitle='time in seconds.', ytitle='no of spikes' ENDELSE ENDELSE end

/ POSSIBLE IMPROVEMENTS: This doesn't have a constant max y-axis value for the graph. So if the max

pro DRAW_CGRAM2, yes_no_values, numberl, number2, binsize, cell_no, time, xmin, xmax this procedure will make a correlogram between the two files input, the functions it uses are :findgen (size) *bin_size+mintemp openfile, find_distance

!p.multi=0

write=string (format= ' ("cross-correlogram of cell",/, (iO) ) ', numberl) write2=string (format= ' (" vs cell ",/,(i0))', number2) timel=where (yes_no_values (numberl -1, *) gt 0) _/this finds where the spikes are time2=where (yes_no_values (number2 -1, *) gt 0) IF ((timel(O) eq -1) OR (time2 (0) eq -1)) then begin mess=WIDGET_MESSAGE ( 'One of the cells has no spikes') ENDIF ELSE BEGIN

N_l=total (yes_no_values , 2) N=N_1 [numberl - 1] ,-this finds the total number of spikes for cell numberl sizeofyn_val=size (yes_no_values) actual__frame_no=sizeofyn_yal [2]

Si=intarr (actual_frame_no) for t=0, actual_frame_no - 1 do begin

Si [t] =yes_no_values [numberl -1 , t] endfor

Sj =intarr (actual_frame_no) for t=0 , actual_frame_no - 1 do begin

Sj [t] =yes_no_values [number2 -l,t] endfor

nlim=fix(- actual_frame_no / binsize) ,-+l _/be careful - some data is lost when tha actual fr.# is not plim=fix(actual_frame_no / binsize) ,--1 ,-be careful - divisible by binsize

,-print, 'nlim - lower bound =', nlim, ' plim - upper bound = ' , plim R=fltarr(2, plim - nlim +1)

index=0 for l=nlim, plim do begin z=0 for t=0, actual_frame_no -1 do begin x=0 for k=t+l*binsize, t+ (1+1) *binsize -1 do begin ,-be careful with the -1 if (k ge 0) and (k It actual_frame_no) then begin x=x+Sj [k] endif endfor if (x gt 0) then begin x=l endif z=z + Si [t] *x _/ (x is either 0 or 1 depending whether cell j fired in a given time bin) endfor R[0, index] =1 R[l, index] =z index=index+1 endfor for i=0, index -2 do begin R[l,i]=R[l,i] / (binsize * N) _/print, R[l,i] endfor

center=-l for i=0, index-2 do begin if (R[0,i] eq 0) then center=i endfor

xminl=fix (xmin/time) xmaxl=fi (xmax/time)

IF binsize It 1 THEN BEGIN mess=WIDGET_MΞSSAGΞ( 'Binsize must be >=1') ENDIF ELSE BEGIN

IF (center + xminl It 0) or (center + xmaxl gt index -1) then begin mess=Widget_Message ( 'X is out of range¹) Endif else begin

xaxis= [ [0, center + xminl] *time, (INDGEN(- xminl + xmaxl) +R[0, center+xminl] +1) *time] yaxis=fltarr (xmaxl-xminl +1) for i=center + xminl, center + xmaxl do begin yaxis [i - center - xminl] = [1, i] endfor plot, xaxis, yaxis , psym=lθ, title=write (0) +write (l)+write2 (0) +write2 (1) , $ xtitle='time in seconds.', ytitle=' firing rate' endelse ENDELSE

ENDELSE end WIDGET FOR CROSS

CORRΞLOGRAMS pro draw_cross common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame_no, file_name,$ total_frame_no , time_resolution, x_size, y_size, box_size common markov2 , yes_no_values3 , yes_no_temp _/added by pvc3 base=WIDGET_BASE (/column, title='Cross Correlogram') draw=WIDGET_DRAW(base, xsize=400, ysize=400) numberl=CW_FIELD (base, title= ' Enter Cell Number 1', VALUE=1, UVALUE=2 , $ /RΞTURN_EVENTS, /INTEGER) number2=CW_FIELD(base, title= ' Enter Cell Number 2', VALUE=2 , UVALUE=3 , $

/RETURN_EVENTS, /INTEGER) binsize=CW_FIΞLD(base, title= ' Enter the binsize', VALUE=1, UVALUE=4 , $

/ ETURN_EVΞNTS , /INTEGER) xmin=CW_FIELD(base, title= 'xl=Lower bound' , VALUE=-200, /RETURN_EVENTS , /INTEGER) xmax=CW_FIELD (base, title= 'x2=Upper bound', VALUE=200, /RETURN_EVENTS , /INTEGER) text=WIDGET_TEXT(base, VALUE=string (' otal number of cells you chose : ' , cell_no) )

WIDGΞT_CONTROL, base, /realize holder={numberl:numberl, number2 : number2 , binsize :binsize, draw: draw, xmax:xmax, xmin:xmin}

WIDGET_CONTROL, WIDGΞT_INFO (base, /Child), SET_UVALUE=holder xmanager, ' draw_cross ' , base end

EVENT HANDLER FOR ABOVE

WIDGET pro draw_cross_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file name, $ tot l_frame_no, time_resolution, x_size, y_size, box_size common markov2, yes_no_values3 _/ added by pvc3 stateholder5=WIDGET_INFO (event . handler, /Child) WIDGET_CONTROL, stateholderδ , GET_UVALUE=state

WIDGET_CONTROL, state. draw, GET_VALUE=win_id2 WIDGET_CONTROL, state . numberl , GET_VALUE=number_l WIDGΞT_CONTROL, state . number2 , GΞT_VALUE=number_2 WIDGET_CONTROL, state. binsize, GΞT_VALUE=bin_size WIDGET_CONTROL, state.xmax, GET_VALUE=xmax WIDGET CONTROL, state.xmin, GET VALUE=xmin

IF ( (number_l gt cell_no) OR (number_l It 1) ) OR ( (number_2 gt cell_no) OR (number_2 It 1) ) THEN BEGIN mess=WIDGET_MESSAGE (string (' Invalid cell number entered. You only chose' , byte (cell_no) , ' cells'), /ERROR) ENDIF ELSE BEGIN wset, win_id2 _/changed by pvc3 draw_cgram, yes_no_values3 , number_l, number_2, bin_size, cell_no, time_resolution, $ xmin, xmax ENDELSE end

WIDGET FOR CROSS

CORRΞLOGRAMS pro draw_cross2 common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common markov2, yes_no_values3, yes_no_temp _/added by pvc3 base=WIDGET_BASE (/column, title='Cell firing rate Cross Correlogram¹ ) draw=WIDGET_DRAW(base, xsize=400, ysize=400) numberl=CW_FIELD (base, title= 'Enter Cell Number I¹, VALUE=1, UVALUE=2 , $

/RETURN_EVENTS , /INTEGER) number2=CW_FIELD(base, title= 'Enter Cell Number 2', VALUE=2 , UVALUE=3 , $

/RETURN_EVENTS , /INTEGER) binsize=CW_FIELD(base, title=' Enter the binsize', VALUE=1, UVALUE=4 , $

/RETURN_EVENTS , /INTEGER) xmin=CW_FIELD(base, title='xl=Lower bound' , VALUE=-200, /RETURN_EVENTS , /INTEGER) xmax=CW_FIELD(base, title= 'x2=Upper bound' , VALUE=200, /RETURN_EVENTS , /INTEGER) text=WIDGET_TEXT(base, VALUE=string (' Total number of cells you chose:', cell_no) )

WIDGET_CONTROL, base, /realize holder= (numberl : numberl, number2 : number2 , binsize :binsize, draw: draw, xmax: max, xmin:xmin}

WIDGET_CONTROL, WIDGET_INFO (base, /Child), SET_UVALUΞ=holder xmanager, ' draw_cross2 ' , base end

_/ EVENT HANDLER FOR ABOVE

WIDGET pro draw_cross2_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common markov2, yes_no_values3 ,- added by pvc3 stateholder5=WIDGET_INFO (event .handler, /Child) WIDGET_CONTROL, stateholderδ , GET_UVALUE=state WIDGET_jCONTROL, state. draw, GET_VALUE=win_id2

WIDGET_CONTROL, state .numberl, GET_VALUE=number_l WIDGET CONTROL, state .number2 , GET VALUE=number 2 WIDGET_CONTROL, state .binsize, GET_VALUE=bin_size WIDGET_CONTROL, state.xmax, GET_VALUE=xmax WIDGET_CONTROL, state.xmin, GET_VALUE=xmin IF ( (number_l gt cell_no) OR (number_l It 1) ) OR ( (number_2 gt cell_no) OR (number_2 It 1) ) THEN BEGIN mess=WIDGET_MΞSSAGE (string ( 'Invalid cell number entered. You only chose ', byte (cell_no) , ' cells'), /ERROR) ENDIF ELSE BEGIN wset, win_id2 _/changed by pvc3 draw_cgram2, yes_no_values3 , number_l, number_2 , bin_size, cell_no, time_resolution, $ xmi , xmax ENDELSE end

,- Name .- draw_raster, vikram, its spelled rastEr, Ted, I know, vikram is spelled with a capital V! ,- Synopsis: rasterjplot

_/ Description: This will draw spikes for those points that pass the test.

_/ Working: It sees if the pixel vs time graph has any peaks and then simply draws a line _/ from the x-axis to represent an 'on' or spike.

It draws the plots one on top of another in a raster form.

Modified: 2/16/97 by VSK

Modified by VSK on 7-11-97 to fix bug where last spike was not printed (changed a=bytarr (frame) to a=bytarr (frame+1) ; Modified by Rob just to get the name straight... let's be consistent, eh?

_/ Modified by Rob, 3.29.99 for formatting. pro draw_raster common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size common choose_cell_convert_delta, filepathl a=bytarr (frame_no+l) loadct, 27

!p.multi=[0, 1, cell_no] window, 1, xsize=600, ysize=600, title=strcompress (strmid(systime () , 3, 8))

_/ Hack to print out raster values to file filename=pickfile (/write, file= 'Location_Spikes_#_' ) if (filename eq ' ' ) then begin

MESS=WIDGET_MESSAGE( 'Data will not be saved!', /INFORMATION) cd, ' ' filename= ' rasterdata . txt ' endif openw, outunit, filename, /GET_LUN printf, outunit, FORMAT= ' ("Number of Frames: ",I0," Number of Cells: ", 10) ' ,frame_no,cell_no printf, outunit, ' ' printf, outunit, 'Cell No No Spikes Spike Times (frame) ' for cell_number = 0, cell_no - 1 do begin if (cell_number le 8) then begin _{/ /} 8 because we add one to the cell number when we print printf, outunit, F0RMAT= '("Cell ",I0,":",I," " , $) ' , cell_number+l, total (yes_no_values (cell_number, *) ) endif else begin printf, outunit, FORMAT= '("Cell ",10,":",!," " , $) ' , cell_number+l, total (yes_no_values (cell_number, *) ) endelse ; / code to write out only those frames that have spikes if (total (yes_no_values (cell_number, *) ) ne 0) then begin _/this only draws cells that have activity plot, a, xstyle=4, ystyle=4, ymargin= [0 , 0] , xmargin= [15, 0] _/xstyle and ystyle avoid display of x and y axis respec. see 2-8R xyouts, - (frame_no/10.000) , 0, cell_number+l, charsize=0.8 _/this prints the cell numbers in front of the cells for frames=0, frame_no-l do begin if (yes_no_values (cell_number, frames) ge 1) then begin _/ to read yes_no plots, [frames, frames], [.75, 0], color=3 printf, outunit, FORMAT= '($,1)', frames endif endfor endif ;; code to write yes_no_values instead- i.e. writes the whole damn binary spike train for each cell ,-printf , outunit, ' ' ,-printf, outunit, transpose (yes_no_values [cell_number, *] ) printf, outunit, ' ' ,- newline endfor

CLOSE, outunit FREE_LUN, outunit !p.multi=0 end

Copyright (c) , 1997 by VSK NAME : draw_significance_raster .pro

SYNOPSIS: draw_significance_raster, yes_no_significance

DESCRIPTION: This program looks at the many_one and two_many significance tests to see / whether many cells fire in one frame or a pair of cells fire together in many frames . ,- It then colors the two cases differently on the raster plot. _/ This says nothing about the significance of the correlation. For that, you must look at the

_/ distribution curve or the data file that pops up after calling many_one or two_many ;

_/ MODIFICATION HISTORY:

_/ Created by Vikram Sheel Kumar on June 20, 1997

_; Last Modified by VSK on jul 12, 1997 to fix the spike printing pro draw_significance_raster, yes_no_significance common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, fi1e_name, $ total_frame_no, time_resolution, x_size, y_size, box_size common with_create_dist, cells_active ,- cells_active is the number of active cells in the movie

DRAW THE BASE LINES FOR THE CELLS AND

PLOT THE SPIKES ON THEM a = bytarr (frame no + 1) _/plotting 'a' will make the line on which the spikes will be for each cell loadct, 27

!p.multi=[0, 1, cells_active] ,-this plots cells_active many single lines one above another window, /free, xsize=600, ysize=600, title= 'Significance Raster Plot' _/the title is the date for cell_number = 0, cell_no - 1 do begin if (total (yes_no_significance (cell_number, *) ) ne 0) then begin ,-this only prints cells with activity plot, a, xstyle=4, ystyle=4, ymargin= [0 , 0] , xmargin= [15, 0] ,-xstyle and ystyle avoid display of x and y axis respec. see 2-8R xyouts, - (frame_no / 10.000), 0, cell_number + 1, color=2 ,-this prints the cell numbers in front of the cells for frames = 0 , frame_no - 1 do begin

_/ COLOR THE SPIKES RED IF THEY HAPPEN

TO FULFILL THE CRITERIA if (yes_no_significance (cell_number, frames) ge 1) then begin _/to read yes_no plots, [frames, frames], [.75, 0], color= (yes no_significance (cell_number, frames) -l)*200+2 _/this makes the correlated spikes red and others grey endif endfor endif endfor

PRINTING A DOTTED BLUE LINE OVER THE MANY

CELLS THAT FIRE TOGETHER IN ONE FRAME for frames = 0, frame no - 1 do begin if total (yes_no_significance (*, frames) gt 1) gt 0 then begin plots, [frames, frames], [200,0], color=3 , linestyle=l endif endfor

!p.multi=0 end _/ Name : draw_spikes

_/ Synopsis: draw_spikes (array)

_/ Description: This will draw spikes for those points that pass the test.

; Working: It sees if the pixel vs time graph has any peaks and then simply draws a line from the x-axis to represent an On' or spike. Created by VSK Last Modified on 12/2/96 pro draw_spikes , number_l common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size common share_w_pixel_vs_time_widget_com, y_min, y_max, spike_min, spike_max, x_max, x_min

!p.multi=0 plot, findgen(frame_no) *time_resolution, pixel_array (number_l-l, *) , yrange= [y_min, y_max] , $ xrange= [x_min, x_max] , title=number_l, ytitle='dF/F' , xtitle='time in seconds', FONT=-l for frames=l, frame_no-l do begin _/this doesn't consider a spike at the 1st frame if (yes_no_values (number_l-l, frames) ne 0) then begin plots, [frames*time_resolution, frames*time_resolution] , [spike_tnax, spike_min] endif endfor end

Description: This program plots the deltaF/F values of a single cell per time. pro draw_spikes_many_cell_delta, delta_f_values, cell number, frame_number, time_resolution

!p.multi=0 window, 2, title= ' Superimposed Plots' plot, findgen(frame_number) *time_resolution, delta_f_values (0, *) , yrange= [min (delta_f_values) ,max (delta_f_values) ] , $ title='All Cells', ytitle= 'dF/F ' , xtitle='time in seconds ' xyouts, -6, 0, 1 if (cell_number gt 1) then begin for cell_numb=l, cell_number-l do begin oplot , f indgen ( frame_number) *time_resolution, (delta_f_values (cell_numb, * ) ) a= (delta_f_values (cell_numb , * ) eq max ( (delta_f_values (cell_numb , * ) ) ) ) a=where (a eq 1) xyouts , a , max ( (delta_f_values (cell_numb , * ) ) ) , cell_numb+l endfor endif end

Name : draw_spikes_many_pages

Modified by Vikram Sheel Kumar on Nov 19, 1996

Modified by Vikram Sheel Kumar on June 23, 1997 to make it only print cells with spikes ; Last Modified by VSK on July 9, 1997 to fix the working in the scenario of no active cells pro draw_spikes_many_pages, pixel_array, yes_no_values, cell_no, frame_number, time_resolution common share_w_pixel_vs_time_widget_com, y_min, y_max, spike_min, spike_max

sum_array=intarr (cell_no) _/find which cells are active for counter=0, cell_no-l do begin sum_array (counter) =total (yes_no_values (counter, *) ) endfor active_cells=total (sum_array gt 0) ,-this has the location of the cells that are active if active_cells eq 0 then begin mess=WIDGET_MESSAGE ( 'There are no active cells!', /error) endif else begin active_cell_subscripts=where (sum_array gt 0) _/this is equal to the number of active cells number = (active_cells / 9) + 1 _/number of windows that have to be opened for t = 0 , number - 1 do begin _/open the windows window, t, xsize=500, ysize=500 endfor !p.multi=[0,3,3] _/it'll print three rows of three cells for index=0, active_cells-l do begin wset, (index) / 9 cell_number=active_cell_subscripts (index) plot, findgen(frame_number) *time_resolution, pixel_array (cell_number, *) , yrange= [y_min, y_max] , $ title=cell_number+l, ytitle= 'dF/F' , xtitle='time in seconds ' for frames=0, frame_number-l do begin if (yes_no_values (cell_number, frames) ge 1) then begin plots, [frames*time_resolution, frames*time_resolution] , [spike_max, spike_min] endif endfor endfor endelse

!p.multi=0 /this returns plotting to systems default end ,- Name : filter_median_subtractive

_/ Description: this was Vikram's 'medsub.pro' file, changed slightly to use with MultiCell

7

_/ Parameters: takes an input array (the array of cells x frames dF/F values) , a window size for

,- filtering, and an output array (an empty variable which will be the filtered input passed back) .

7

_/ History: Created February 25, 1999, by Rob from code by Vikram. _/ I thought so! We shouldn't truncate at 0! Changed 3.25.99. pro filter_median_subtractive, input_array, median_window_size, output_array output_array = float (input_array) num_cells = n_elements (input_array (* , 0 ) ) num_frames = n_elements (input_array (0, *) ) temp_array = fltarr (num_frames + (2 * median_window_size) ) for i = 0, (num__cells - 1) do begin temp_array (median_window_size :median_window_size + num_frames - 1) = reform (input_array (i, *) ) temp_array = (temp_array - median (temp_array, median_window_size) ) output_array (i , *) = temp_array (median_window_size :median_window_size + num_frames - 1) endfor end

^~~; Copyright (cl New Ron Productions

_/ Name: find_distance

_/ Synopsis: find_distance, timearrayl, timearray2 _/ Description: This procedure is used by draw_cgram to make cross correlograms and auto correlograms .

_/ Timearrayl and timearray2 are arrays which hold the times where the spikes occur. The output

,- of this procedure is an array of the difference of the second spike train from the first cell ' s .

_/ It does this by^' taking the spike of the first cell eg 13 54 66 154 and looks at the second spike

,- train's timearray which may look like: 14 55 67 155

_/ outputarray (0, 0) =14-13=1 / outputarray ( 0 , l)=14-54=-40

,- outputarray (0, 2)=14-66=-52

_/ outputarray (0, 3) =14-154=-140 outputarray (1, 0) =55-13=42 outputarray (1, 1) =55-54 =1 and so on this will show a peak at 1 in the cross correlogram for 1 to 2.

function find_distance, timearrayl, timearray2

_/this function takes two files and finds how their elements differ by subtracting all of the elements ,- of one from the other. Output is a file with a measure of their deviation from each other sizel=n_elements (timearrayl) size2=n_elements (timearray2) outputarray=intarr (sizel, size2) index=0 for i=0, sizel-1 do begin ,-one by one subtracting the elements of filel

^■ _/ from the whole of file2 outputarray (index, *) =timearray2 -timearrayl (i) index=temporary (index) +1 endfor _/for j=0, sizel-1 do begin

_/for t=0, size2-l do begin _/outputarray (j , t) =timarray2 (t) -timearrayl (j ) _/ endfor _/ endfor return, outputarray end

_/ NAME:

_/ find_matches .pro

7 ,- SYNOPSIS:

_/ find_matches , celll, cell2, winsize

7

_/ INTRODUCTION :

; This calculates the correlation coefficient of two cells with a binsize

7

; INPUTS:

; celll and cell2 are the unidimensional arrays of the timearrays of the cells. _/ winsize is the size of the window. THIS MUST BE AN ODD NUMBER.

_/ DESCRIPTION:

_/ The c coeff is calculated as follows. _/ 1. The binary spikes are converted into time arrays where the values of the

_/ array correspond to the indices where the spikes occurred.

_/ For example, if celll=[0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1,

1] ,- timel=where (celll gt 0)

_/ IDL> print, where (celll gt 0)

4 7 10 11 2. The IDL library function strpos is used to find whether any of the time array values are the same . If there are two time array values that are the same, that means that both of the cells we are analyzing had spikes at the same time. That counts as a match in the formula to calculate correlation coefficients : c_coeff= number of matches (A- >B) / total spikes of A

3. find_matches counts the number of times there are congruent spike positions. The number of these occurences = matches

4. The total number of spikes in cell 1 =lengthl

5 Cross correlation coef=matches/lengthl

MODIFICATION HISTORY

Created May 1996 by Vikram Sheel Kumar

Last Modified Jul 13, 1997 to remove a redundant step

MAIN PROGRAM

BEGINS function find_matches, celll, cell2, winsize matches=0 ,- matches is the number of coincident spikes timel=where (celll gt 0) ; gives the location of where celll is greater than 0 time2=where (cell2 gt 0) gt means greater than, similarly- >ge is greater than equal lengthl=total (ti el gt 0) lengthl gives the number of elements in the array timel

;which is also the number of spikes in celll the same window (whose size is defined by winsize) as a spike from cell2. It is calculated in the module below

FIND THE NUMBER OF COINCIDENT

SPIKES for i=0, lengthl-1 do begin window_size= (winsize-1) /2 _/a window is created x before a spike and x after. Therefore making it 2x+l large lower_window_limit=timel (i) -window_size upper_window_limit=timel (i) +window_size count=0 _/count is the number of times a spike from celll occurs in for m=lower_window_limit, upper_window_limit do begin temp=strpos (time2, m) templ= (temp) +1 _/there is a zero if there is a match therefore we have to add 1 if (total (tempi) ne 0 ) then count=count+l _/ in this case there is a coincident spike endfor if (count gt 0) then matches=matches+l endfor

_/ CALCULATE THE CORRELATION

COEFFICIENT coeff=matches/ (lengthl) ,- lengthl is always greater than 0 as we only do this for active cells return, [coeff] end

gen_sig_widget

This is the widget for the General Significance test.

Created, Rob, 4.21.99

Gen sig test revised to include windowing and remove many/once and two/many, 1.00, Rob .

pro gen_sig_widget gen_sig_base=WIDGΞT_BASE (/COLUMN, title= ' General Significance') gen_sig_num_trials=CW_FIELD (gen_sig_base , /RETURN_ΞVENTS ,

/INTEGER, TITLE= ' Number of random trials to run:', VALUE=1000,

UVALUE=0) gen_sig_window_size=CW_FIELD(gen_sig__base, /RETURN_EVENTS,

/INTEGER, TITLE='Window size:', VALUE=0, UVALUE=0) gen_sig_button=WIDGΞT_BUTTON (gen_sig_base , VALUE= ' Find significance', UVALUE=2) WIDGET_CONTROL, /realize, gen_sig_base gen_sig_state={gen_sig_num_trials :gen_sig_num_trials, $ gen_sig_window_size :gen_sig_window_size, $ gen_sig_button : gen_sig_button} WIDGET_CONTROL, WIDGET_INFO (gen_sig_base, /CHILD), SET_UVALUE=gen_sig_state xmanager, 'gen_sig_widget ' , gen_sig_base end

1997, Copyright by Vikram Sheel Kumar Name: gen_sig_widget_event Synopsis: gen_sig_widget_event

Description: this program sits in Multicell and calls the other statistical tests. It calls

_/ find_general_p which gives a general statistical measure of the analyzed movie. It also looks at

_/ the analyzed movie to see how many cells fire in one frame and how many _/ time two cells fire together in the movie. _/ Modified by Rob, 4.10.99, to add all the code from

' find_general_p ' here just because it seemed useless to have it as two separate files.

_/ Modified, Rob, 4.21.99, turning this into the event handler like the Many Once and Two Many tests

_/ Modified, Rob, 4.24.99, speedup for Monte Carlo- ignores silent cells

_{/ /} Gen sig test revised to include windowing and remove many/one and two/many, 1.00, Rob. pro gen_sig_widget_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with_create_dist, cells_active gen_sig_stateholder=WIDGET_INFO (event .handler, /CHILD) WIDGΞT_CONTROL, gen_sig_stateholder, GET_UVALUΞ=gen_sig_state WIDGET_C0NTR0L, gen_sig_state . gen_sig_num_trials ,

GET_VALUE=no_of_times

WIDGET_C0NTR0L, gen_sig_state .gen_sig_window_size,

GET_VALUE=window_size if (event. id eq gen_sig_state.gen_sig_button) then begin if ( (no_of_times le 1) or (window_size le -1) ) then begin error_dialog = WIDGETJV1ESSAGE ( ' Invalid fields specified! ' , /error) endif else begin / ; We write all of our data into the following file in the function find_general_p filename3 = 'General_Statistical_Data' _/this is what Multicell xdisplays, so don't change name filename3 = pickfile (/write, file='General_Stats_#') if (filename3 eq ' ' ) then begin mess = WIDGET_MΞSSAGE ( ' The Data will not be saved.', /INFORMATION) filename3 = 'General_Statistical_Data' /this is what Multicell xdisplays, so don't change name endif

WIDGET CONTROL, /HOURGLASS

/ ,- creating the random distribution times_repeated = 2 no_spikes = total (yes_no_values) distribution_array = intarr (no_of_times) _{/ /} array of hits for each random trial for t = 0 , no_of_times - 1 do begin make_random_data, seed, random_array, num_active_cells distribution_array (t) = count_random_match.es (random_array, times_repeated, window_size) endfor _{/ /} count the number of matches in the real data now true_match.es = 0 if (window_size eq 0) then begin for frame_counter = 0, frame_no - 1 do begin matches = total (yes_no_values (*, frame_counter) ) /this gives us the n in our combination formula above if (matches ge times_repeated) then begin matches_pair = factorial (matches)

/ (factorial (times_repeated) * factorial (matches - times_repeated) ) ,- finds number of possible pairs endif else begin matches_pair = 0 ,-cant find factorial of a negative number endelse true_matches = true_match.es + matches_pair /increases the matches count for the true data by that of this frame endfor endif else begin

/ ,- first compress spike trains according to window size binned_yn_vals_array = intarr (cell_no, frame_no - (window_size * 2)) for frame_counter = window_size, frame_no - (window_size + 1) do begin for window_counter = -window_size, window_size do begin for cell_counter = 0, cell no - 1 do begin if (yes_no_values [cell_counter, frame_counter + window_counter] ge 1) then begin binned_yn_vals_array [cell_counter, frame_counter - window_size] = 1 endif endfor endfor endfor

/ ,- then look for matches in the binned spike trains (although a single spike can fall into many bins) for frame_counter = 0, n_elements (binned_yn_vals_array [0, *] ) - 1 do begin matches = total (binned_yn_vals_array (*, frame_counter) ) /finds number of hits at this frame if (matches ge times_repeated) then begin matches_pair = factorial (matches)

/ (factorial (times_repeated) * factorial (matches - times__repeated) ) ,- finds number of possible pairs endif else begin matches_pair = 0 ,-cant find factorial of a negative number endelse true_match.es = true_match.es + matches_pair _/ increases the matches count for the true data by that of this frame endfor endelse

,- ,- determine connectivity stats distances_array = [0.0] total_connections = 0 connections_array = intarr (cell_no) if (window_size eq 0) then begin for cell_l = 0, cell_no - 2 do begin for cell_2 = cell_l + 1, cell_no - 1 do begin temp_cell = yes_no_values (cell_l,

*) * yes_no_values (cell_2, *) if ( (total (temp_cell) ) ge 1) then begin connections_array [cell_l] = connections_array [cell_l] + 1 connections_array [cell_2] = connections_array [cell_2] + 1 total_connections = total_connections + 1 ,- ,- do this so a connection is only counted once, if it matters connection_distance = sqrt ( ( (double (location (cell_l) . coord [0] ) - double (location (cell_2) .coord [0] ) ) 2) + ( (double (location (cell_l) . coord [1] ) - double (location (cell_2) .coord [1] ) ) 2)) distances_array = [distances_array, connection_distance] endif endfor endfor endif else begin for cell_l = 0 , cell_no - 2 do begin for cell_2 = cell_l + 1, cell_no - 1 do begin temp_cell = binned_yn_vals_array (cell_l, *) * binned_yn_vals_array (cell_2 , *) if ( (total (temp_cell) ) ge 1) then begin connections_array [cell_l] = connections_array [cell_l] + 1 connections_array [cell_2] = connections_array [cell_2] + 1 total_connections = total_connections + 1 _/ ; do this so a connection is only counted once, if it matters connection_distance = sqrt (( (double (location (cell_l) . coord [0] ) - double (location (cell_2) .coord [0] ) ) 2) + ( (double (location (cell_l) . coord [1] ) - double (location (cell_2) .coord [1] ) ) ^Λ 2)) distances_array = [distances_array, connection_distance] endif endfor endfor endelse

_{/ /} Begin printing the data close, 1 openw, 1, filename3 printf, 1, 'STATISTICAL DATA' printf, 1, • ' printf, 1, 'Random Test Information¹ printf, 1, ' ' printf, 1, 'How many times a pair is together: ', strcompress (times_repeated) printf, 1, 'Number of iterations in tests:

' , strcompress (no_of_times) printf, 1, 'Window size: ' , strcompress (window_size) printf, 1, 'Total number of active cells: ', strcompress (round (cells_active) ) printf, 1, ' ' if (max(distribution_array) eq 0) then begin mess = WIDGET_MESSAGE ( ' Two cells never fire together randomly! Try the MANY CELLS AT ONCE option looking at 2 cells firing together with more than 100 iterations.', /error) printf, 1, 'ERROR! ' printf, 1, 'Two cells never fire together randomly. Nothing can be said about the significance of this data' real_p = -1 endif else begin loadct, 27 hist = histogram (distribution_array, binsize=l, min=0, max= (ma ( [true_matches,distribution_array] ) + 2)) window, /free, title= 'Distribution for

General Significance' plot, hist, xtitle= 'number of hits', ytitle= ' frequency' plots, [true_matches, true_matches] , [0, max ( [1, total (distribution_array eq true_matches) ] ) ] , color=12

_/this draws a line where the actual number of matches lies spikes_per_cell_per_second = (no_spikes / cells active) / (frame_no * time resolution) stats = moment (distribution_array, sdev=gen_sdev) p_value = gauss_pdf ( (true_matches - stats (0))

/ gen_sdev) real_p = 1 - p_value printf, 1, 'General Statistical Information' printf, 1, ' ' printf, 1, 'Total number of frames: ', strcompress (frame_no) printf, 1, 'Time between frames: strcompress (time_resolution) printf, 1, 'Total number of cells: strcompress (cell_no) printf, 1, 'Total number of spikes: strcompress (no_spikes) printf, 1, 'Mean expected matches: strcompress (stats (0) ) printf, 1, 'Variance: strcompress (stats (1) ) printf, 1, 'Standard deviation: strcompress (gen_sdev) printf, 1, 'Actual matches: ', strcompress (true_match.es) printf, 1, 'Actual/expected: ' , strcompress (true_matches / stats (0)) printf, 1, 'Standard error for ratio:

', strcompress (gen_sdev / stats (0)) printf, 1, 'Mean spikes per cell: ', strcompress (no_spikes / cells_active) printf, 1, 'Spike firing rate: ' , strcompress (spikes_per_cell_per_second) printf, 1, 'Significance p-value: ', strcompress (real_p) printf, 1, ' (IMPORTANT! If this is exactly .5, it could be actually very significant so look at actual/expected!) ' printf, 1, ' ' printf, 1, 'Cell Connectivity Information' printf, 1, ' ' printf, 1, 'Total Number of Connections:

' , strcompress (total_connections) connected_cells = where (connections_array) num_connected_cells = n_elements (connected_cells) if (n_elements (connected_cells) ge 2) then begin stats_conn_only = moment (connections_array (connected_cells) , sdev=sdey_conn_only) mean_connections = stats_conn_onl [0] normalized_slice_connectivity = (mean_connections / num_connected_cells) printf, 1, 'Mean Number of Connections per Cell: ', strcompress (mean_connections) printf, 1, 'Standard Deviation:

¹ , strcompress (sdev_conn_only) endif else begin normalized_slice_connectivity = 0.0 ,- ,- no connected cells! Can't have just one connected cell. endelse printf, 1, 'Normalized Connectivity: ' , strcompress (normalized_slice_connectivity) printf, 1, 'Number of Silent Cells: ' , strcompress (cell_no - num_connected_cells) printf, 1, 'Number of Coactive Cells:

' , strcompress (n_elements (connected_cells) ) printf, 1, 'Minimum Number of Connections: strcompress (min(connections_array) ) printf, 1, 'Maximum Number of Connections: strcompress (max(connections_array) ) if (n_elements (distances_array) gt 2) then begin distance_stats = moment (distances_array [1 :n_elements (distances_array) - 1], sdev=sdev_distances_array) printf, 1, 'Mean Connection Distance: ' , strcompress (distance_stats [0] ) printf, 1, 'Standard Deviation: ' , strcompress (sdev_distances_array) endif endelse

_{/ /} Close the file we are writing, reopen it if small enough close, 1 openr, 1, filename3 stat_stuff = fstat(l) file_size = stat__stuff . size close, 1 if (file_size gt 8112) then begin mess=WIDGET_MESSAGE( 'File is too large to display through a widget. Open it manually. If you did not save it, it is named General_Statistical_Data' , /INFORMATION) endif else begin xdisplayfile, filename3 , title = "General

Statistical Data", group = event. top, width = 75, height = 50 endelse endelse endif end load_and_convert_buqing. ro

This is for loading in and converting one of Bu Qing's data files.

Created 5.08.99, Rob.

5.13.99, reads in diode positions and builds maps

Warning! Currently, this depends on 60x magnification (and hence tries to fit a 12.3 micron distance between diodes) . _{/ /} Also, the number of diodes, 464, is hardcoded. pro load_and_convert_textfile common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common test, str_ing, state3 common old_skool_data, original_data common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined common with_choose_cells_com, thresholdl, threshold2 , rms_threshold common textfile_vars, text_flag, diode_array, max_num_of_diodes max_num_of_diodes = 464 message_dialog=WIDGET_MESSAGE ("First, please select a file to read the individual diode traces from. An IDL ascii template will then come up, please first select the line the data begins at, press 'Next', then type the number of diodes, press 'Next', then just press ' Finish ' . " , /INFORMATION) filename=dialog_pickfile (/read, file= (' 990401fl-l.txt ') , get_path=filepathl) if (filename eq ' ' ) then begin message_dialog=WIDGET_MESSAGE("No data read.", /INFORMATION) endif else begin close, /all openr, 1, filename

_{/ /} opens one of IDL's stupid ascii templates so that the user can select where the data begins

_{/ /} and how many frames there are. Why isn't FRAME a keyword to 'read_ascii ' ??

_{/ /} anyways, 'header' returns the first four lines- two of which are null strings, the other two which we ; _/ use to determine the time resolution and the diode numbers a little further down in the code... data = read_ascii (filename, template=ascii_template (filename) , data_start=4, header=header) close, /all cell_no = n_tags (data) keyword_pos = strpos (header [0] , 'Number of Data Points Per Trace: ' ) frame no = fix(strmid (header [0] , keyword_pos + 32, 6)) - 1

_//frame no = n_elements (data. fieldOl)

; _/ throw all the crap away to get the time resolution and the photodiode numbers keyword_pos = strpos (header [0] , 'Frame Interval:') time_resolution = double (strmid (header [0] , keyword_pos + 15)) / 1000.0 diode_array = intarr (cell_no) current_position = 0 for i = 0, cell_no - 1 do begin det_pos = strpos (header [2] , 'Det') comma_pos = strpos (header [2] , ' , ' ) diode_array [i] = fix(strmid (header [2] , 3 + det_pos, comma_pos - 1) ) header [2] = strmid (header [2] , comma_pos + 2) endfor

WIDGET_CONTROL, state3.name, SET_VALUE=strmid (filename, strlen (filepathl) ) WIDGET_CONTROL, state3. frame, SET_VALUE=frame_no

WIDGET CONTROL, state3.time, SET_VALUE=time_resolution _/ ,- build 'pixel_array' pixel_array = intarr (cell_no, frame_no) for i = 0, cell_no - 1 do begin pixel_array [i, *] = data. (i) endfor original_data = pixel_array message_dialog=WIDGET_MΞSSAGE ("Now select the 'abscissa' (x-position) file.", /INFORMATION) filename=dialog_pickfile (/read, file= ( ' pda_abscissa . txt ' ) , get_path=filepathl) if (filename eq ' ' ) then begin message_dialog=WIDGΞT_MESSAGE("No data read.", /INFORMATION) endif else begin close, /all openr, 1, filename endelse x_pos_data = read_ascii (filename, data_start=0, header=header)

_/ ,- here we throw out all the extra zeros- we skip the first two data points (a string and a leading 0)

_{/ /} and take the next 'max_num_of_diodes ' datapoints . x_pos_data = x_pos_data. fieldl [2 :max_num_of_diodes + 1] ,- ,- '+ 1' because we do -1 (to read from 0-423) and then +2 (to skip the first two elements in the array) close, /all message_dialog=WIDGET_MΞSSAGE( "Finally, select the 'ordinate' (y-position) file.", /INFORMATION) filename=dialog_pickfile (/read, file= ( ' pda_ordinate . txt ' ) , get_path=filepathl) if (filename eq ' ' ) then begin message_dialog=WIDGET__MESSAGE ("No data read.", /INFORMATION) endif else begin close, /all openr, 1, filename endelse y_pos_data = read_ascii (filename, data_start=0, header=header) y_pos_data = y_pos_data. fieldl [2 :max_num_of_diodes + 1] _/ ,• ditto as above close, /all

_/ / magnification correction : compare 1 , 2 _/ 236 , 237 _/ 100 , 101 _/ 341 , 342 _/ 228 , 229 _/ and 463 , 464

; ; remember ! 12 . 3 is hardcoded for 60x magnification ! We also assume correction is the same in both x and y direction ! sample_set = [x_pos_data [0] - x_pos_data [1] , x_pos_data [235 ] - x_pos_data [236] , x_pos_data [99] - x_pos_data [100 ] , $ y_pos_data [0] - y_pos_data [1] , y_pos_data [235] - y_pos_data [236] , y_pos_data [99] - y_pos_data [100] , $ x_pos_data [340 ] - x_pos_data [341] , x_pos_data [227] - x_pos_data [228] , x_pos_data [462 ] - x_pos_data [463 ] , $ y_pos_data [340] - y_pos_data [341] , y_pos_data [227 ] - y_pos_data [228] , y_pos_data [462 ] - y_pos_data [463 ] ] sample_set_stats = moment (sample_set) correction_factor = 12.3 / sample_set_stats [0] x_pos_data = x_pos_data * correction_factor y_pos_data = y_pos_data * correction_factor

_/ ,- initialize other variables yes_no_values = intar (cell_no, frame_no) coef = fltarr (cell_no, cell_no) rms_threshold = fltarr (cell no) thresholdl = fltarr (cell_no) threshold2 = fltarr (cell_no) for cell = 0, cell no - 1 do begin rms_threshold[cell] = 2.0 thresholdl [cell] = 2.0 threshold [cell] = 3.0 endfor cells_defined = 1 spikes_defined = 0 single_plot_defined = 0 correl_coef_defined = 0 box_size = 2 x_size = max (x_pos_data) + 350 y_size = max(y_pos_data) + 350

,- ,- locations are shifted up and to the right by 300 so that all locations are positive (and thus can be plotted correctly) location = replicate ({struct, coord: intarr (2) , size:0, half_side:0.00} , cell_no) for i = 0, (cell_no - 1) do begin location (i) . size = cell_no location (i) .half_side = 1

_{/ /} error checking for missing diodes- put them right in the middle, why not? if (diode_arra [i] eq 0) then begin location (i) . coord [0] = 300 location (i) . coord [1] = 300 endif else begin location (i) .coord [0] = x_pos_data [diode_array [i] - 1] + 300 location (i) .coord [1] = y_pos_data [diode_array [i] - 1] + 300 endelse endfor text_flag message_dialog=WIDGΞT_MESSAGE( "Finished reading data.",

,- Name : make_binary

,- Synopsis: make_binary, pixel_delta_values, frame_no, cell_no, thresholdl, threshold2

_/ Description: This program will find the spikes after either the pixel intensities of the

_/ cells have been measured or the deltaF values . It is a common function that comes under Threshold. _/ Created by Vikram Sheel Kumar

; Modified on 12/3/96 by VSK to make spikes be made at the increasing or decreasing points (not at the increased or decreased points) ,- Modified on 3.12.97 by VSK to be able to calculate the threshold for an individual cell, (called make_single_binary.pro) _/ Modified, 2.99, by Rob to make threshold vars into arrays function make_binary, delta_values, frame_no, cell_no, thresholdl, threshold2 binary_values = intarr (cell no, frame_no) for cell = 0, cell_no - 1 do begin Case 1 of

(thresholdl [cell] gt 0): begin _/for increasing spikes, the threshold will be greater than 0 for frame=l, frame_no-l do begin ,-it doesn't consider a spike at the first frame if (frame It frame_no-2) then begin

_/ for all of the peaks upto the last frame if ( ( (delta_values (cell, frame+1) ) - (delta_values (cell, frame)) ge thresholdl [cell] ) OR $

(delta_values (cell, frame+2)- delta_values (cell, frame) ge threshold2 [cell] ) ) then begin binary_values (cell, frame) =1 endif endif else begin ,-this takes care of a peak at last frame if ( (delta_values (cell, frame_no-l) ) -

(delta_values (cell, frame_no-2) ) ge thresholdl [cell] ) then begin binary_values (cell, frame_no-2) =1 endif endelse endfor end

(thresholdl [cell] It 0): begin _/for decreasing spikes, the threshold will be less than 0 for frame=l, frame_no-l do begin ,-it doesn't consider a spike at the first frame if (frame It frame_no-2) then begin if ( ( (delta_values (cell, frame+1) ) - (delta_values (cell, frame)) le thresholdl [cell] ) OR $ (delta_values (cell, frame+2) - delta_values (cell, frame) le threshold2 [cell] ) ) then begin binary_values (cell, frame) =1 endif endif else begin _/this takes care of a peak at last frame if ( (delta_values (cell, frame_no-l) ) - (delta_values (cell, frame_no-2) ) le thresholdl [cell] ) then begin binary_values (cell, frame_no-2) =1 endif endelse endfor end endcase endfor return, binary_values end

1997, (c) Copyright by Vikram Sheel Kumar .... Name : make_random_data .pro Synopsis: make_random_data, seed

Description: This program looks at the input cell's activity and creates a random set of data based on this data, by rotating each cell's spike train by a random amount .

Created on May 21, 1997

Modified by Rob, 2.99, to fix the infinite loop bug and the spikes-aren' t-being-shuffled-right ' bug...

4.10.99, Rob, deleted code that was commented out just to keep file small. If you want to see what was here before, look at some of the backup source. 4.24.99, Rob- fundamental change and speedup- this is no longer function, it's a procedure that returns both the number of active cells and a random array of those cells' spike trains, randomly rotated.

Please note that this random array contains ONLY those cells which spike at least once! pro make_random_data, my_seed, random_array, num_active_cells common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with_create_dist, cells_active

/ find the number of spikes each cell of the true data has no_spikes = intarr (cell_no) ,-this is the array with the number of spikes for each cell for cell = 0, cell_no - 1 do begin no_spikes (cell) = total (yes_no_values (cell, *) ) endfor cells_active = total (no_spikes gt 0) num_active_cells = cells_active random_array = intarr (cells_active, frame_no) random_cell_index = 0 for cell_counter = 0, cell_no - 1 do begin if (no_spikes (cell_counter) ge 1) then begin

,- ,- pick a random int between 0 and frame_no random_time_shift = fix(randomu(my_seed) * frame_no) random_array [random_cell_index, *] = shift (yes_no_values [cell_counter, *] , random_time_shift) random_cell_index = random_cell_index + 1 endif endfor end Copyright (c) 1996, Brain Is Ltd. Name: make_single_binary

Synopsis : make_single_binary, delta_values, yes_no_values, cell_number, frame_no, cell_no, threshl, thresh2

Description:

This program comes under the single plot call to Multicell. It is called when the user enters 2 thresholds that

,- are different from the 2 thresholds that calculated the spikes for the whole movie (under the choose threshold call to

,- multicell) . Make_binary calculated the spikes for the whole movie while this function is called make_single_binary.

; This calculates the spikes only when told to do so for a single cell under the single plot widget. Therefore it is called _/ make_single_binary.

_/ It goes through the deltaF/F values (delta_values) and sees if the absolute increase in delta_value between adjacent frames ,- is greater than or equal to the threshold for positive spikes and less than or equal to the threshold for negative spikes.

We do not consider spikes at the first frame.

Explanation of variables : delta_values : the array with the actual delta F/F values yes_no_values : the converted binary array. cell_number: the cell whose spikes are being found. frame_no : total number of frames threshl: threshold between two adjacent frames thresh2 : threshold between three adjacent frames

Modification History

Created by Vikram Sheel Kumar on 3.12.97

Last Modified by VSK on 7/9/97 to exclude finding a spike in the first frame. function make_single_binary, delta_values , yes__no_values, cell number, frame_no, threshl, thresh2 yes_no_values (cell_number-l, *) =0 Case 1 of

(threshl gt 0) : begin _/for increasing spikes, the threshold will be greater than 0 for frame=l, frame_no-l do begin _/this doesn't consider a spike at the first frame if (frame It frame_no-2) then begin if ( ( (delta_values (cell_number-l, frame+1) ) - (delta_values (cell_number-l, frame)) ge threshl) OR $

(delta_values (cell_number-l, frame+2) - delta_values (cell_number-l, frame) ge thresh2) ) then begin yes_no_values (cell_number-l, frame) =1 endif endif else begin _/this takes care of a peak at point 1 if ( (delta_values (cell_number-l, frame_no-l) ) - (delta_values (cell_number-l, frame no-2)) ge threshl) then begin yes_no_values (cell_number-l, frame_no-2) =1 endif endelse endfor end

(threshl It 0): begin _/for decreasing spikes, the threshold will be less than 0 for frame=l, frame_no-l do begin /this doesn't consider a spike at the first frame if (frame It frame_no-2) then begin if ( ( (delta_values (cell_number-l, frame+1) ) - (delta_values (cell_number-l, frame)) le threshl) OR $

(delta__values (cell_number-l, frame+2) - delta_values (cell_number-l, frame) le thresh2)) then begin yes_no_values (cell_number-l, frame) =1 endif endif else begin _/this takes care of a peak at point 1 if ( (delta_values (cell_number-l, frame_no-l) ) -

(delta_values (cell_number-l, frame_no-2) ) le threshl) then begin yes_no_values (cell_number-l, frame_no-2) =1 endif endelse endfor end endcase return, yes_no_values end many cells_one_frame

Creates a random distribution (via count_random_matches_many and make_random_data)

_{/ /} to compare the real data to- counts the number of times a minimum number of

; ,- cells fire in one frame, or a set of frames given by the window size.

77

,- _/ Originally written, Vikram, 6.2.97 _{/ /} Rewritten, Rob, 5.26.99 (with miscellanious improvements in- between) function many_cells_one_frame, synchronous_cells, no_iterations, ge or eq test, window_size common mother_com, pixel_array, yes_no_values, coef, location, cell_no , frame_no, fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size common with_create_dist, cells_active common with riany, filename

_/common with_correl_map_plan, cluster_array ,-uncomment this if you want networks to have different linestyles in the correl_map_plane.. _{/ /} Initialize variables connections_array = bytarr (cell_no, cell_no) _/ if two cells are significantly connected, this is 1_/ if not, 0 yes_no_significance = yes_no_values _/ the new colored raster plot is made from yes_no_significance . This sets it to the original data cluster_array = intarr (cell_no, cell no) _/ this holds the frame in which the cells fire together true_hits = 0 _/ number of networks in the real data distances_array = [0.0] _/ list of all network connection distances areas_array = [0.0] _/ list of all network areas size_array = intarr (frame_no) _/ this counts only those networks ' synchronous_cells '+ in size cluster_no = 0 _; this makes the first network in solid lines

_{/ /} Creating the random distribution no_spikes = total (yes_no_values) _/this give the total number of spikes by cells in the slice random_distribution = intarr (no_iterations) _/this array holds one random number of matches for every iteration for t = 0, no_iterations - 1 do begin make_random_data, seed, random_array, num_active_cells random_distribution(t) = count_random_matches_many (random_array, synchronous_cells, ge_or_eq_test , window_size) endfor

_/ ,- Printing some of the data close, 1 filename = pickfile (/write, file= (strcompress (string (synchronous_cells) ) + '_Cells_Once ' ) ) if (filename eq ' ' ) then begin mess = WIDGET_MESSAGE ( ' This data will not be saved ' , /INFORMATION) filename = 'Many_One.dat' endif openw, 1, filename _/DO THIS-this is the temp file used for writing into the xdisplay. It can be saved to a further name . printf, 1, 'Statistical data for the Many/Once test.' printf, 1, 'One hit is when ', strcompress (synchronous_cells) , ' cells fire in one frame' printf, 1, 'Number of iterations: ', no_iterations printf, 1, 'Window size: ', window_size printf, 1, 'Below is the list of cells which fire in a given frame : '

,- ,- Finding the real number of matches- first bins spikes according to window size, no binning done when window is 0 binned_yn_vals_array = intarr (cell_no, frame_no - (window_size * 2) ) for frame_counter = window_size, frame_no - (window_size + l) do begin for window_counter = -window_size, window_size do begin for cell_counter = 0, cell no - 1 do begin if (yes_no_values [cell_counter, frame_counter + window_counter] ge 1) then begin binned_yn_vals_array [cell_counter, frame_counter - window_size] = 1 endif endfor endfor endfor _{/ /} Now looks for matches in the binned data- there may be false positives for windows for frame_counter = 0, n_elements (binned_yn_vals_array [0, *] ) - 1 do begin matches = total (binned_yn_vals_array (* , frame_counter) ) if ( ( (ge_or_eq_test eq 0) and (matches ge synchronous_cells) ) or ( (ge_or_eq_test eq 1) and (matches eq synchronous_cells) ) ) then begin true_hits = true_hits + 1 size_array [frame_counter] = matches printf, 1, ' ' if (window_size eq 0) then begin printf, 1, strcompress (round (matches) ) , ' Cells Fire in Frame: ', strcompress (frame_counter) correlated_cell_array = where (yes_no_values (*, frame_counter) eq 1) endif else begin printf, 1, strcompress (round (matches) ) , ' Cells Fire at Approximately Frame: ', strcompress (frame_counter) correlated_cell_array = where (binned_yn_vals_array (*, frame_counter) eq 1) endelse network_xcors = intarr (n_elements (correlated_cell_array) ) network_ycors = intarr (n_elements (correlated_cell_array) ) for index = 0, matches - 1 do begin printf, 1, ' Cell number: ', correlated_cell_array (index) + 1 loc = correlated_cell_array (index) for window_counter = frame_counter, frame_counter + (2 * window_size) do begin if (yes_no_values (correlated_cell_array(index) , window_counter) ge 1) then begin yes_no_significance (correlated_cell_array (index) , window_counter) = yes_no_significance (correlated_cell_array (index) , window_counter) + 1 endif endfor network_xcors [index] = location (loc) . coord [0] network_ycors [index] = location (loc) .coord [1] printf, 1, ' X and Y: ' , location (loc) .coord printf, 1, ' ' for index2 = index, matches - 1 do begin if (index ne index2) then begin connection_distance = sqrt ( ( (double (location (correlated_cell_array [index] ) . coord [0] ) - double (location (correlated_cell_array [index2] ) .coord [0] ) ) 2) + ( (double (location (correlated_cell_array [index] ) . coord [1] ) - double (location (correlated_cell_array [index2] ) .coord [1] ) ) ^Λ 2)) distances_array = [distances_array, connection_distance] endif loc2 = correlated_cell_array (index2) connections_array (loc, loc2) = 1 cluster_array (loc, loc2) = cluster_no endfor endfor

_{/ /} find area of many cells once network if ( (n_elements (uniq (network_xcors , sort (network_xcors) ) ) ge 3) and (n_elements (uniq(network_ycors, sort (network_ycors) ) ) ge 3) ) then begin triangulate, network_xcors , network_ycors , triangles , convex_hull / _/ here we find the convex hull surrounding the network convex_hull = [convex_hull, convex_hull [0] ] _/ the last vertex is the first for connections ' sake area = 0.0D for i = 0, (n_elements (convex_hull) - 2) do begin area = area +

( (double (network_xcors [convex_hull [i] ] ) * double (network_ycors [convex_hull [i + 1] ] ) ) - (double (network_xcors [convex tmll [i + 1] ] ) * double (network_ycors [convex_hull [i] ] ) ) ) endfor area = abs (area / 2) areas_array = [areas_array, area] printf, 1, ' Area of network (in pixels^Λ2) : ', strcompress (area) printf, 1, ' ' endif cluster_no = cluster no + 1 endif endfor

_{/ /} draw raster plot before it checks to see if the random distribution has a variance as the raster is independant of statistics draw_significance_raster, yes_no_significance

,- ,- draw correlations map window, /free, xsize=x_size, ysize=y_size, title='Correl Map of All Many/Once Networks' _/ connections_array has the data of the cells connected correl_map_image_plane, connections_array

,- _/ draw histogram of areas if (n_elements (areas_array) ge 2) then begin window, /free, title= 'Distribution of Many/Once Areas' hist_areas = histogram (areas_array [l:n_elements (areas_array) - 1], binsize=5000 , min=0, max=max(areas_array) + 5000) plot, lindgen (n_elements (hist_areas) + 2) * 5000L, hist_areas, psym=10, title= 'Distribution of Many/Once Areas ', $ yrange=[0, max(hist_areas) + 1], xrange=[0, max(areas_array) + 5000], ytitle= 'Number of Networks', xtitle='Area (Pixels 2) ' endif

/ / check to see if the random distribution has any non zero values. If not, the moment cannot be defined for the distribution xmax=max (random_distribution) xmin=min (random_distribution) if (xmax eq xmin) then begin _/if the min and max are the same, all of the elements are the same mess=WIDGET_MΞSSAGE( 'Random distribution has variance of zero. Try again with greater number of iterations', /error) printf, 1, 'ERROR!' printf, 1, 'Moment undefined for random distribution with variance zero ' close, 1 endif else begin

_/plot random distribution with a line for the actual value

Ip.multi = 0 bin_size = 1 hist = histogram (random_distribution, binsize = 1, min =

0, max ! * xmax + 1) ,-plot a histogram of the distribution window, /free, title='Many One Distribution' plot, hist, xtitle= 'number of hits', ytitle= 'frequency ' y2 = total (random_distribution eq true_hits) if y2 le 0 then begin ,-if there is no random data=true_hits, to draw the blue line y2 = y2 + 1 endif plots, [true_hits, true_hits] , [0, y2] , color=12 ,-this draws a line where the actual number of matches lies

_/calculating the number of spikes per cell per second spikes_per_cell_per_second = (no_spikes / cells_active) / (frame no * time resolution)

_/calculating the p value, standard dev. for the data stats = moment (random_distribution, sdev=sdev) no_points_right = total (random_distribution ge true_hits) p_value = no_points_right / no_iterations

_/writing more data to the file: printf, 1, 'Total number of active cells: strcompress (cells_active) printf, 1, 'Total number of frames: strcompress (frame_no) printf, 1, 'Total number of spikes: ', strcompress (no_spikes) printf, 1, 'Mean expected matches: ', strcompress (stats (0) ) printf, 1, 'Variance: ' strcompress (stats (1) ) printf, 1, 'Standard deviation: ' strcompress (sdev) printf, 1, 'Actual matches (no. of networks): ', strcompress (true_hits) printf, 1, 'Number of networks/number of cells: ', strcompress (true_hits / cells_active) printf, 1, 'Normalized number of networks: ' strcompress (true_hits / cells_active / frame_no) printf, 1, 'Actual/expected: ', strcompress (true_hits / stats (0)) printf, 1, 'Standard error for ratio: ' strcompress (sdev / stats (0)) if (n_elements (where (size_array) ) ge 2) then begin network_stats = moment ( size_array (where (size_array) ) , sdev=sdev_size_array) printf, 1, 'Mean cells in a network: ', strcompress (network_stats [0] ) printf, 1, 'Standard deviation: ', strcompress (sdev_size_array) endif printf, 1, ' ' printf, 1, 'Spike firing rate: ' strcompress (spikes_per_cell_per_second) if (n_elements (distances_array) gt 2) then begin distance_stats = moment (distances_array [1 :n_elements (distances_array) - 1], sdev=sdev_distances_array) printf, 1, 'Mean connection distance: ' , strcompress (distance_stats [0] ) printf, 1, 'Standard deviation: ' , strcompress (sdey_distances_array) endif else begin if (n_elements (distances_array) gt 1) then begin printf, 1, 'Connection distance:

' , strcompress (distances_array [1] ) endif endelse if (n_elements (areas_array) gt 2) then begin area_stats = moment (areas_array [1 :n_elements (areas_array) - 1], sdev=sdev_areas_array) printf, 1, 'Mean area: ' strcompress (area_stats [0] ) printf, 1, 'Standard deviation:

' , strcompress (sdey_areas_array) endif printf, 1, 'Significance p-value: ', strcompress (p_value) close, 1 endelse return, [random_array] end

(c) , 1997 Copyright by Vikram Sheel Kumar

Name : many_one_widget

Synopsis : many_one_widget

Description: Through this widget, we can find how significant many cells firing in one frame is. The user

_/ can choose how many cells consititutes a hit and how many iterations the program should run through to create the random_distribution.

Modified on 4.20.99, Rob, to add radio buttons for 'eq' or 'ge' testing control _/ Modified, 5.8.99, Rob, to add window code pro many_one_widget common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size many_one_base=WIDGET_BASE (/COLUMN, title=' Significance of many cells firing together in one frame') many_one_fieldl=CW_FIELD (many_one_base , /RETURN_ΞVENTS , /INTEGER, TITLE= ' Number of cells to fire in one frame to count as a hit:',$

VALUE=4, UVALUΞ=0) ge_or_eq_bgroup=CW_BGROUP (many_one_base, ['This Many or More' , 'Exactly This Many'] , /ROW, /EXCLUSIVE, SET_VALUE=0) many_one_field2=CW_FIELD(many_one_base, /RETURN_EVENTS ,

/INTEGER, TITLE= 'Number of iterations :', $

VALUΞ=1000, UVALUE=0) window_size_field=CW_FIELD (many_one_base , /RETURN_ΞVΞNTS , /INTEGER, TITLE= ' Window size for hits:', VALUE=0, UVALUE=0) many_one_buttonl=WIDGET_BUTTON (many_one_base, VALUE= ' Find significance¹, UVALUE=2)

WIDGET__CONTROL, /realize, many_one_base many_one_state={many_one_fieldl :many_one_fieldl, $ many_one_field2 :many_one_field2 , $ window_size_field:window_size_field, $ many_one_button :many_one_buttonl , $ ge_or_eq_bgroup : ge_or_eq_bgroup} WIDGET_CONTROL, WIDGΞT_INFO (many_one_base, /CHILD), SΞT_UVALUE=many_one_state xmanager, ' many_one_widget ' , many_one_base end

©Copyright Synapse Innovation Inc .

Written by Vikram Sheel Kumar ,- NAME: many_one_widget_event , event

SYNOPSIS: many_one_widget_event , event (called by xmanager)

DESCRIPTION: This is the event handler for the many_one_widge .

Created on June 18, 1997 by Vikram Sheel Kumar

Modified on July 26, 1997 by VSK to check for size of file before xdisplaying

_/ Modified on 4.20.99, Rob, to add radio buttons for 'eq' or 'ge' testing control _/ Modified, 5.8.99, Rob, to add window code pro many_one_widget_event , event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with_many, filename stateholder3=WIDGET_INFO (event .handler, /CHILD) WIDGET_CONTROL, stateholder3 , GET_UVALUE=many_one_state

WIDGET_CONTROL, many_one_state .many_one_fieldl , GET_VALUE=synchronous_cells

WIDGET_CONTROL, many_one_state .many_one_field2 , GET_VALUE=no_of_iterations WIDGΞT_CONTROL, many_one_state . indow_size_field, GET_VALUE=window_size

WIDGET_CONTROL, many_one_state . ge_or_eq_bgroup , GET_VALUE=ge_or_eq_test if (event. id eq many_one_state.many_one_button) then begin if ( (synchronous_cells le 0) OR (no_of_iterations le 1) ) then begin mess=WIDGET_MΞSSAGE( 'Invalid fields specified!', /error) , endif else begin

I random_array = many_cells_one_frame (synchronous_cells , no_of_iterations , ge_or_eq_test, window_size)

WIDGET_CONTROL, event . top , /hourglass close, 1 openr, 1, filename stat_stuff = fstat(l) file_size = stat_stuff .size close, 1 if (file_size gt 8112) then begin mess = WIDGET_MΞSSAGΞ ( ' File is too large to display through a widget. Open it manually. If you did not save it, it is named Many_One.dat', /INFORMATION) endif else begin xdisplayfile, filename, title = "Statistical

Data for Many/One", group = event. top, width = 75, height = 50 endelse endelse endif end Created, VSK

,- Modified, Rob Froemke, 2.99... changed default frame number because it was pissin' me off,

2.22.99, added wave correlation interface ; 3.16.99, added save/load variables to/from textfiles

3.17.99, init all common block vars ... 5.13.99, stripped down interface, removing extraneous widgets and controls. ,- 1.00, revised gen sig test pro MultiCell common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common old_skool_data, original_data common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined common with_choose_cells_com, thresholdl, threshold2, rms_threshold common test, str_ing, state3 common choose_cell_convert_delta, filepathl common textfile_vars, text_flag, diode_array, max_num_of_diodes

_{/ /} Init all the common block vars above! (except state3) pixel_array = [0] yes_no_values = [0] original_data = [0] coef = [0] location = 0 cell_no = 1 frame_no = 1 f ile_name = ' ' total_frame_no = 300 time_resolution = 1.0 x_size = 0 y_size = 0 box_size = 0 cells_defined = 0 spikes_defined = 0 single_plot_defined = 0 correl_coef_defined = 0 rms_threshold = [0] thresholdl = [0] threshold2 = [0] str_ing = ' Movie . ' buqing_flag= 0 diode_array = [0] filepathl = 'dev/null' _/ ,- but where else to send our data? text_flag = 0 diode_array = [0]

base=WIDGET_BASE (/column, title= 'Multicell ' )

_/ ,- setting uvalue=0 below tells the event handler that we're going to refer to the buttons by name,

_/ ,- i.e. pass the button name as a string as the event . b=cw_bgroup (base, /row, [ ' ew Experiment ' , $ ' Save as XDR ' , $ ' Load from XDR ' , $ 'Exit'], /return_name, UVALUE=0)

_/this is the structure thai z hhas the information for the pulldown menu stuff= {cw_pdmenu_s , flags :0, name : ' ' } details= [ { cw_pdmenu_s, 1, ' Find Spikes ' } , $

{ cw_pdmenu_s , o, ' Root Mean Squared ' } , $

{ cw_pdmenu_s , 2, ' Intensity Threshold ' } , $

{ cw_pdmenu_s, 0, 'Delete Spikes ' } , $

{ cw_pdmenu_s , 1, 'Plotting' }, $ cw_pdmenu_s , 0 , 'Single Plots' } , $ cw_pdmenu_s , 0 , 'All Plots' },$ cw_pdmenu_s , 0 , 'All PlotS-3D' },$ cw_pdmenu_s , 0 , 'Superimposed Plots '}, $ cw_pdmenu_s , 0 , 'Raster Plot' },$ cw_pdmenu_s , 2 , 'Overall Behavior'}, $ cw_pdmenu_s , 0 , 'Load from Text File'}] pull_down=cw_pdmenu(base, details, /return_full_name, UVALUE=12) stuff2= {cw_pdmenu_s2, flags :0, name:''} details= [ cw_pdmenu_s2 , 1 'Test Significance'}, $ cw_pdmenu_s2 , 0 'General Significance'}, $ cw_pdmenu_s2 , 0 'Many Cells One Time'}, $ cw_pdmenu_s2 , 2 'Two Cells Many Times'}, $ cw_pdmenu_s2 , 0 'Build Correlation Map ' } , $ cw_pdmenu_s2, 0 ' Cross Correlogram ' } ] pull_down2=cw_pdmenu(base, details, /return_full_name, UVALUE=13) stuff3= {cw_pdmenu_s3, flags :0, name:''} details= [ { cw_pdmenu_s3, 0, 'Color Tables¹}, $ { cw_pdmenu_s3 , 0, 'Return to IDL'}] pull_down3=cw_pdmenu (base, details, /return_full_name, UVALUE=14) namesid=cw_field(base, title= ' Filename Root ', UVALUE=2, VALUE=str_ing, /STRING, /return_events) frameid=CW_FIELD(base, title= ' otal Number of Frames ' ,VALUE=total_frame_no, /INTEGER, /RETURN_EVENTS) timeid=CW_FIELD (base, title='Time Resolution (seconds) ' ,VALUE=time_resolution, /FLOATING, /RETURN_EVENTS) widget_control , /realize, base ,- make the widget visible state={name amesid, frame :frameid, time:timeid} WIDGET_CONTROL, WIDGET_INFO (base, /Child), SET_UVALUE=State xmanager, 'MultiCell', base /register this widget with the xmanager so it can call the event end ,- Name: MultiCell_event

,- Synopsis: MultiCell_event, event

_/ Description: This is the event handler for the MultiCell widget.

It registers all events with

,- the xmanager and calls the various procedures that need to be called.

_/ Created by Vikram Sheel Kumar on June 23, 1997 as an adaptation of

Mother4.0

,- Modified on July 26, 1997 by VSK to check for size of file before xdisplaying a file under general. stats... _/ Modified, 2.99, Rob, threshold arrays: we need to do a back- compatibility check when we load variables...

_/ And again, to add the ' correlate_all_waveforms ' control. Modified 3.16.99, buttons for textfile saving/loading 5.13.99, complete redesign- extraneous controls eliminated. 1.00, revised gen sig test

pro MultiCell_event, event

,- ,- important... if you want your variables saved in the textfile, make sure to add code to ' save_to_text_file ' and ;; ' load_from_text_file¹ ! common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common old_skool_data, original_data common test, str_ing, state3 common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined common with_choose_cells_com, thresholdl, threshold2 , rms_threshold common choose_cell_convert_delta, filepathl common markov2 , yes_no_values3 , yes_no_temp , cell_no_temp , frame_no_temp common textfile_vars, text_flag, diode_array, max_num_of_diodes stateholder3=WIDGET_INFO (event. handler, /Child) WIDGET_CONTROL, stateholder3 , GET_UVALUΞ=state3 ,-find what is in /Child

WIDGET_CONTROL, state3.name, GET_VALUE=file_name WIDGET_CONTROL, state3. frame, GET_VALUE=total_frame_no WIDGET CONTROL, state3.time, GET VALUE=time resolution

if ( (size (event .value) ) [1] eq 7) then begin _// i.e. if we've pressed a button which returns a string as its value... case event.value of

-FIRST ROW ITEMS-

¹ New Experiment ' : begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (total_frame_no le 4) then begin mess=WIDGET_MESSAGE('I need at least 5 frames to calculate the baseline!', /INFORMATION) endif else begin choose_cells_chief endelse endelse end

'Save as XDR': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin data_file=pickfile (/write, file= ' Variables_from_exp_#_' , title= ' Create the saved variables file') if (data_file eq ' ' ) then begin mess=WIDGET_MΞSSAGE ( 'No saved variables file specified.', /INFORMATION) endif else begin

SAVE, /VARIABLES, FILENAME=data_file, all, /verbose mess=WIDGET_MESSAGE('Data from this experiment saved!', /INFORMATION) endelse endelse end

'Load from XDR' : begin state4=state3 if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MΞSSAGE( 'Enter frame number and time resolution before proceeding! ' , /ERROR) endif else begin data_file=pickfile (/read, title= ' Select the saved variables file', GET_PATH=filepathl) if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGE ( 'No saved variables file specified. ' , /INFORMATION) endif else begin rms_threshold = [0] _{/ /} reset it so that it can be properly reinitialized if need be restore, data_file size_of_loc=size (location) temp_cell=size_of_loc (1) size_of_y_n_v=size (yes_no_values) frame_no=size_of_y_n_v (2 ) cell__no=size_of_y_n_v (1) if (temp_cell ne cell no) then begin mess=WIDGET_MESSAGE (' Something is wrong with your saved variables file!', /Error) endif else begin

/ / cells_defined = 1 ; / spikes_defines = 1 single_plot_defined=0 text_flag = 0 if (n_elements (thresholdl) eq 1) then begin _{/ /} back-compatibility check tempi = thresholdl temp2 = threshold2 thresholdl = fltarr (cell_no) threshold2 = fltarr (cell_no) for cell = 0 , cell_no

- 1 do begin thresholdl [cell] = tempi threshold2 [cell] = temp2 endfor endif _/ ; more back- compatibility... if

(n_elements (rms_threshold) eq 1) then begin rms_threshold = fltarr (cell_no) for cell = 0, cell_no

- 1 do begin rms_threshold [cell] = 2 . 0 endfor endif

WIDGET_CONTROL, state4. name , SET_VALUE=str_ing

WIDGET_CONTROL, state4. frame , SΞT_VALUE=frame_no

WIDGET_C0NTR0L, state4. time , SET_VALUE=time_resolution mess=WIDGET_MESSAGΞ ( ' Data from previous experiment loaded!', /INFORMATION) endelse endelse data_file=0 endelse end

'Exit' :begin

WIDGET_CONTROL, /DESTROY, event . top end

-SECOND ROW ITEMS-

'Find Spikes. Root Mean Squared': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (cells_defined ne 1) then begin mess=WIDGET_MESSAGE('You have to find cells before you find spikes ! ' , /error) endif else begin rms_spikes_widget endelse endelse end

'Find Spikes . Intensity Threshold': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (cells_defined ne 1) then begin mesε=WIDGET_MESSAGE ( 'You have to find cells before you find spikes ! ' , /error) endif else begin choose_threshold endelse endelse end 'Delete Spikes': begin if (spikes_defined ne 1) then begin mess = WIDGET_MESSAGE ( 'You have to find spikes before you can delete any! ' , /error) endif else begin delete_spikes_widget endelse end

' Plotting. Single Plots': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceedin ! ' , /ERROR) endif else begin if (spikes_defined ne 1) then begin mess=WIDGET_MESSAGE ( 'You have to find spikes before you can plot cells ! ' , /Error) endif else begin pixel_vs_time_widget endelse endelse end

'Plotting.All Plots': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE( 'Enter frame number and time resolution before proceeding! ' , /ERROR) endif else begin if (single_plot_defined ne 1) then begin mess=WIDGET_MESSAGΞ('You have to plot single cells before you can plot all cells ! ' , /Error) endif else begin draw_spikes_many_pages , pixel_array, yes_no_values, cell no, frame_no, time_resolution endelse endelse end 'Plotting.All Plots-3D': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (cells_de ined ne 1) then begin mess=WIDGET_MESSAGE( 'You have to plot single cells before you can plot all cells ! ' , /Error) endif else begin draw_3D_plot, pixel_array endelse endelse end

'Plotting. Superimposed Plots': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (cells_defined ne 1) then begin mess=WIDGET_MESSAGE(' You must find cells before you can plot them! ' , /Error) endif else begin draw_spikes_many_cell_delta, pixel_array, cell_no, frame_no, time_resolution endelse endelse end

'Plotting. Raster Plot¹: begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MΞSSAGE( 'Enter frame number and time resolution before proceeding!¹, /ERROR) endif else begin if (spikes_defined ne 1) then begin mess=WIDGET_MESSAGE('You have to find spikes before drawing a raster plot ! ' , /Error) endif else begin if (total (yes_no_values) eq 0) then begin mess=WIDGET_MESSAGE ( 'None of the cells have any spikes! Try lowering the thresholds', /INFORMATION) endif else begin draw_raster endelse endelse endelse end

'Plotting. Overall Behavior': begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE( 'Enter frame number and time resolution before proceeding!¹, /ERROR) endif else begin if ( single_j?lot_def ined ne 1) then begin mess=WIDGET_MESSAGΞ('You have to plot single cells before you can plot all cells ! ' , /Error) endif else begin summed_spikes endelse endelse end

'Load from Text File' : begin load_and_convert_textfile end / / THIRD ROW ITEMS

'Test Significance. General Significance' : begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess = WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding! ' , /ERROR) endif else begin if (spikes_defined ne 1) then begin mess = WIDGET_MESSAGE ( 'You have to find spikes before calculating the significance of correlations ! ' , /error) endif else begin gen_sig_widget endelse endelse end

'Test Significance.Many Cells One Time':begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess = WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (spikes_defined ne 1) then begin mess = WIDGET_MESSAGE ( 'You have to find spikes before calculating the significance of correlations ! ' , /error) endif else begin many_one_widget endelse endelse end

'Test Significance. Two Cells Many Times' :begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess = WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (spikes_defined ne 1) then begin mess = WIDGET_MESSAGE ( 'You have to find spikes before calculating the significance of correlations ! ' , /error) endif else begin two_many_widget endelse endelse end 'Build Correlation Map':begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess = WIDGET_MESSAGE (' Enter frame number and time resolution before proceeding ! ' , /ERROR) endif else begin if (spikes_defined ne 1) then begin mess = WIDGET_MESSAGE ( 'You have to find spikes before calculating correlation coefficients ! ' , /error) endif else begin widget_analyze endelse endelse end 'Cross Correlogram' :begin if (total_frame_no eq 0) or (time_resolution eq 0) then begin mess=WIDGET_MESSAGE ( 'Enter frame number and time resolution before proceeding!', /ERROR) endif else begin if (spikes_defined ne 1) then begin mess=WIDGET_MESSAGE('You have to find spikes before you can draw correlograms ! ' , /Error) endif else begin yes_no_values3 = yes_no_values choose_crosscorr endelse endelse end

-FOURTH ROW ITEMS-

' Color Tables ': begin xloadct end

'Return to IDL' :begin retail end else: endcase endif end _/Name: pick_box Version4.0 _/Synopsis: picks_box, cell_numbers

_/Description: Through this function, the user picks the cells he/she wants . These coordinates can then

,-be input into the get_cells_oplot_win function which gives the pixel values and binary transformation of the cells.

UPDATED VERSION AS OF 6/25/96 by Vikram Sheel Kumar ,-this program is a version of pick_cells that highlights a cell as soon as a point in the cell is _/clicked.

,-picks_box can only be done after pick_cells has been done so the boxes can be _/ chosen from the same image as the cells were chosen

7 _/OLD statements that might be imp. later: naming the files : -

NOTE, CAN'T READ FILES WITH LAST NAME 'TIFF' names=realname (0) +strcompress ( (strin (format=' (i) ' , findgen(fram e_number) +1) ) , /remove_all)

7 names=filename (0) +strcompress ( (string (format= ' (i) ' , findgen (fram e_number)+l) ) , /remove_all)

_/ names=filename (0) +strcompress ( (string (format= ' (i, " . tiff") ' , findgen (frame_number) +1) ) , /remove_all) _/ data being saved: -

_/ file_name2=string ( ("Saved_Locations_Box@") , strcompress (strmid (systime () , 3 , 10) , $

/remove_all) , '_' , strcompress (strmid (systime () , 14, 2) ) ) close, 3 openw, 3, FILEPATH(file_name2, subdirectory=dir) _/OPEN LUN TO WRITE INTO

; printf, 3, 'half side of box= ' , box_size

_/ printf, 3, 'For cell ', i+1, 'you clicked at {', cell (i) .coord, '}^■

,- close, 3 _/CLOSE THE LUN

MAIN PROGRAM

BEGINS function pick_box, box_size, xsize, ysize

celll=replicate ( {struc, coord: intarr (2) } , 500) _/Stru to hold the location (x and y) of MAX NO OF CELLS CHOOSABLE=500

/ for the first cell boxlength=box_size*2 cell_number=l cursor, x, y, /up, /device celll (cell number) . coord (0)=x celll (cell number) . coord (l)=y

. give an error if the mouse was not clicked on the image if ( (x ge xsize-boxlength OR x le boxlength) OR (y ge ysize- boxlength OR y le boxlength) ) then begin mess=WIDGET_MESSAGE( ' ou clicked outside the allowed area! Type xmanager to return to Multicell.', /Error) retail endif

,- to make sure the first cell doesn't fail the while test we make cell(O) diff from cell(l) celll(O) .coord(0)=-l for every other cell while (x ne celll (cell_number-l) . coord (0) OR y ne celll (cell number-1) .coord (1) ) do begin

,-if these cond fail, you have double clicked on a cell

_/and the program stops counting cells.

. draw the box and number for j=-box_size, box_size do begin plots, [x+j , x+j], [y-box_size, y+box_size] , /device /plots lets me draw points endfor xyouts, x-boxlength*3, y+boxlength, cell_number, $ ,-xyouts lets you write on an image alignments 5, /device cell number=cell number+1 read the cursor position from the screen cursor, x, y, /up, /device give an error if the mouse was not clicked on the image- if ( (x ge xsize-boxlength OR x le boxlength) OR (y ge ysize- boxlength OR y le boxlength) ) then begin mess=WIDGET_MESSAGE( 'You clicked outside the allowed area! Type xmanager to return to Multicell.', /Error) retail endif

/ put the x and y coordinates into the structure celll (cell_number) .coord (0) =x celll (cell_number) . coord (l)=y endwhile make another structure with the correct no of cells i.e cell number-1 cell=replicate ( {struct, coord : intarr (2) , size:0, half_side : 0.00} , cell_number-l) for i=0, cell_number-2 do begin cell (i) . coord=celll (i+1) .coord endfor cell . size=cell_number-l _/total number of cells is cell_number-l since one extra cell is cell.half_side=box_size /recorded when you double click to stop counting. return, [cell] end

,- while not doing so for the whole array of cell .

,- This will work by having the threshold widget work as it did before, and then also _/ creating an option to rethreshold particular cells if necessary.

,- Through this method, the spikes will be recalculated every time the graph is displayed.

_/ Or it can have a statement which checks whether the threshold chosen for that cell is the ,- universal one, and in that case it will not rethreshold.

; Created on 3.12.97 by VSK to be able to calculate indiv. thresholds for different cells.

_/ Modified, 2.99, by Rob for lots of stuff (threshold arrays, spike/cell deletion, initial plotting, 1 cell, ; automatic y-ranging, location of textboxes)

_/ 3.25.99, Rob, delete trailing spikes in a spike clump

_/ 5.9.99, Rob, complete redesign for real estate and efficiency, also adding individual cell RMS thresholds,

,- individual cell smoothing/filtering/restoration _/ WIDGET FOR DRAWING A SINGLE

SPIKE pro pixel_vs_time_widget common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with__choose_cells_com, thresholdl, threshold2, rms_threshold deep_base=WIDGET_BASE(/row, title='Plot of a Single Cell', /scroll) left_base=WIDGET_BASE(deep_base, /column) right_base=WIDGET_BASE (deep_base, /column) draw=WIDGET_DRAW(left_base, xsize=450, ysize=450) ;; kludge to fix the 'only one cell selected' bug if (cell_no gt 1) then begin slide=WIDGET_SLIDER(left_base, value=l, maximum=cell_no, minimum=l, $ title= ' Choose the cell you want to plot', uvalue= ' slider_pressed' ) endif else begin slide=WIDGET_TEXT(left_base, value= ' Only one cell selected... ' ) endelse spike_max=CW_FIELD(right_base, title='Y Max for Spike:', VALUE=0, /RΞTURN_EVΞNTS , /FLOATING) spike_min=CW_FIELD (right_base, title='Y Min for Spike: ', VALUE=-20, /RETURN_EVENTS , /FLOATING) ymax=CW_FIELD(right_base, title= 'Y Range Max:', VALUE=max( [transpose (pixel_array (0, *) ) , 10]), /RETURN_EVENTS, /FLOATING) ymin=CW_FIELD (right_base, title='Y Range Min: ', VALUE=min( [transpose (pixel_array ( 0 , *) ) , -10]), /RETURN_EVENTS , /FLOATING) xmax=CW_FIELD (right_base, title='X Range Max:', VALUE=total_frame_no*time_resolution, / ETURN_EVENTS , /FLOATING) xmin=CW_FIΞLD(right_base, title='X Range Min: ', VALUE=0, /RETURN_EVENTS, /FLOATING) median_filter_window=CW_FIELD (right_base, title= 'Window for Median Filter: ', VALUE=20, /RETURN_EVENTS , /FLOATING) smoothing_window=CW_FIELD (right_base, title= 'Window for Mean Smoothing:', VALUE=3 , /RETURN_ΞVENTS , /FLOATING) buttonsl=cw_bgroup (right_base, /row, ['Smooth', 'Median Filter' , 'Restore Original Waveform'], $

BUTTON_UVALUE= [ ' smooth_pressed' , ' filter_pressed' , ¹ restore_pressed' ] , UVALUE=5) rms_threshold_box=CW_FIELD(right_base, title='RMS Threshold:', VALUE=2.0 , $

UVALUE=19, /RETURN_EVENTS , /FLOATING) rms_spikes_type_bgroup=CW_BGROUP (right_base , [ ' Positive Spikes' , 'Negative Spikes'] , /ROW, /EXCLUSIVE, SET_VALUE=0) threshold_one=CW_FIELD (right_base, title= ' 2 Frame Intensity Diff Threshold:', VALUE=thresholdl [0] , $

UVALUE=6, /RETURN_EVENTS , /FLOATING) threshold_two=CW_FIELD (right_base, title= ' 3 Frame Intensity Diff Threshold:', VALUE=threshold2 [0] , $

UVALUE=7, /RETURN_EVENTS , /FLOATING) buttons2=cw_bgroup (right_base, /row, ['Plot with RMS', 'Plot with Intensity Difference'], $

BUTTON_UVALUΞ= [ ' RMS_plot_pressed ' , ' diff_plot_pressed ' ] , UVALUE=3 ) delete_spike=CW_FIELD (right_base, title= 'Delete Spike Number:¹, VALUE= ' 0 ' , UVALUE=2 , $

/RETURN_EVENTS , /INTEGER) buttons3=cw_bgroup (right_base, /row, ['Delete Cell', 'Delete All Trailer Spikes ' ] , $

BUTTON_UVALUE= [ ' delete_cell_pressed' , 'delete_trailer_spikes_pressed' ] , UVALUE=4) holder={draw:draw, slide:slide, delete_spike :delete_spike, $ spike_max : spike_max, spike_min : spike_min, $ ymax:ymax, ymin:ymin, xmin:xmin, xmax:xmax,$ rms_spikes_type_bgroup : rms_spikes_type_bgroup , $ median_filter_window:median_filter_window, smoothing_window : smoothing_window, $ rms_threshold_box : rms_threshold_box, threshold_one : threshold_one , threshold_two : threshold_two} WIDGET_CONTROL , deep_base , /realize WIDGET ONTROL, deep_base , set_uvalue=holder plot, findgen (frame_no) * time_resolution, pixel_array (0, *) , yrange= [min( [transpose (pixel_array(0, *)), -10]), ma ( [transpose (pixel_array (0, *) ) , 10])],$ xrange=[0, total_frame_no * time_resolution] , title=l, ytitle='dF/F' , xtitle='time in seconds', FONT=-l for frames=l, frame_no - 1 do begin _/this doesn't consider a spike at the 1st frame if (yes no_values (0, frames) ne 0) then begin plots, [frames * time_resolution, frames * time_resolution] , [0, -20] endif endfor

XMANAGER, 'pixel_vs_time_widget ' , deep_base end

,- Created by Vikram S . Kaze

; Modified by Rob Froemke, 2.99- lots of stuff (threshold arrays, spike/cell deletion, cell-by-cell yranges)

_/ 3.25.99: control for trailing spike clump deletion ; 5.9.99: more redesigning... individual RMS thresholds pro pixel_vs_time_widget_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common share_w_pixel_vs_time_widget_com, y_min, y_max, spike_min, spike_max, x_max, x_min common with_choose_cells_com, thresholdl, threshold2, rms_threshold common old_skool_data, original_data common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined WIDGET_CONTROL, event. top, get_uvalue=state

WIDGET_CONTROL, event. id, GET_UVALUE=uval

WIDGET_CONTROL, state. draw, GET_VALUE=window

WIDGET_CONTROL, state . slide , GET_VALUE=number_l if (string (uval) eq ' sliderjpressed' ) then begin

WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0

WIDGET_C0NTR0L, state.ymax, set_value = max (pixel_array ( (number__l - 1), *) ) y_max = ma (pixel_array ( (number_l - 1), *) )

WIDGET_CONTROL, state.ymin, set_value = min (pixel_array ( (number__l - 1), *) ) y_min = min (pixel_arra ( (number_l - 1), *) ) _{/ /} check for 10,-10 boundary... cheap...

WIDGET_CONTROL, state.ymax, set_value = max([y_max, 10]) y_max = max([y_max, 10]) WIDGET_CONTROL, state.ymin, set_value = min([y_min, -10]) y_min = min([y_min, -10])

WIDGET_CONTROL, state . rms_threshold_box, SET_VALUE=rms_threshold[number_l - 1] rms_t = rms_threshold [number_l - 1] WIDGET_CONTROL, state . threshold_one, SET_VALUE=thresholdl [number_l - 1] tl = thresholdl [number_l - 1] WIDGET_CONTROL, state . threshold_two,

SET_VALUE=threshold2 [number_l - 1] t2 = threshold2 [number_l - 1] endif WIDGET_CONTROL, state . delete_spike, GET_VALUE=spike_l

WIDGET_CONTROL, state .median_filter_window, GET_VALUE=median_filter_window

WIDGET_CONTRO , state . smoothing_window, GET_VALUE=smoothing_window WIDGET_CONTROL, state . rms_threshold_box, GET_VALUE=rms_t

WIDGET_CONTROL, state . threshold_one , GET_VALUE=tl

WIDGET_CONTROL, state . threshold_two, GΞT_VALUE=t2

WIDGET_CONTROL, state.ymin, GET_VALUE=y_min

WIDGET_CONTROL, state.ymax, GET_VALUE=y_max WIDGET_CONTROL, state.xmin, GET_VALUE=x_min

WIDGET_CONTROL, state. max, GET_VALUE=x_max

WIDGΞT_CONTROL, state . spike_min, GET_VALUE=spike_min

WIDGET_CONTROL, state . spike_max, GΞT_VALUE=spike_max wset, window _/even if you open other windows, this will still plot in the original window

,- _/ trailing spike clumps deletion if (string (event .value) eq 'delete_trailer_spikes_pressed' ) then begin for frame_counter = (frame_no - 1), 2, -1 do begin if ( (yes_no_values [number_l - 1, frame_counter] eq 1) and ( (yes_no_values [number_l - 1, frame_counter - 1] eq 1) or (yes_no_values [number_l - 1, frame_counter - 2] eq 1))) then begin yes_no_values [number_l - 1, frame_counter] =

0 endif endfor draw_spikes, number_l endif

_// if we're re-plotting, update the yes_no_values to current threshold if (string (event .value) eq ' diff_plot_pressed' ) then begin WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0 if ( (tl eq 0) OR (t2 eq 0)) then begin mess=WIDGET_MESSAGE( 'Threshold cannot be 0!', /Error) WIDGET_CONTROL, state . threshold_one,

SET VALUE=thresholdl [number_l - 1] WIDGET CONTROL, state . threshold_two, SET_VALUE=threshold2 [number_l - 1] tl = thresholdl [number_l - 1] t2 = threshold2 [number_l - 1] endif else begin thresholdl [number_l - 1] = tl threshold2 [number_l - 1] = t2 yes_no_values=make_single_binary (pixel_array, yes_no_values, number_l, frame_no, thresholdl [number_l - 1], threshold2 [number_l - 1] ) endelse endif if (string (event .value) eq 'RMS_plot_pressed' ) then begin WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0

WIDGET_CONTROL, state . rms_spikes_type_bgroup, GET_VALUE=pos_or_neg rms_threshold [number_l - 1] = rms_t _{/ /} update new threshold with value in textbox if (pos_or_neg eq 0) then begin cell_stats = moment (pixel_array [number_l - 1,*], sdev=cell_dff_stddev) cell_dff_mean = cell_stats [0] for frame_counter = 1, frame_no - 1 do begin if (pixel_array [number_l - 1, frame_counter] gt (cell_dff_mean + (rms_threshold [number_l - 1] * cell_dff_stddev) ) then begin yes_no_values [number_l - 1, frame_counter] = 1 endif else begin yes_no_values [number_l - 1, frame_counter] = 0 endelse endfor ,- ,- negative spikes is pos_or_neg = 1 endif else begin cell_stats = moment (pixel_array [number_l - 1,*], sdev=cell_dff_stddev) cell_dff_mean = cell_stats [0] for frame_counter = 1, frame_no - 1 do begin if (pixel_array [number_l - 1, frame_counter] It (cell_dff_mean - (rms_threshold[number_l - 1] * cell_dff_stddev) ) then begin yes_no_values [number_l - 1, frame_counter] = 1 endif else begin yes_no_values [number_l - 1, frame_counter] = 0 endelse endfor endelse endif

_{/ /} filtering- subtractive median filter: see filter_median subtractive. pro if (string (event.value) eq ' filterjoressed' ) then begin if ( (median_filter_window It 2) or (median_filter_window ge frame_no) ) then begin _//error checking error nessage = WIDGET_MESSAGE ( 'Window size for filtering must be greater than 1 and less than the number of frames ! ' , /INFORMATION) endif else begin filter_median_subtractive, pixel_array [number_l - 1,*], median_filter_window, filtered_cell pixel_array [number_l - 1,*] = filtered_cell endelse endif

_// smoothing, via the IDL function 'smooth' if (string (event .value) eq ' smooth_pressed' ) then begin if ( (smoothing_window It 2) or (smoothing_window ge frame_no) ) then begin _//error checking error_message = WIDGET_MESSAGE ( 'Window size for smoothing must be greater than 1 and less than the number of frames!', /INFORMATION) endif else begin pixel_array [number_l - 1,*] = smooth (pixel_array [number_l - 1,*], smoothing_window) endelse endif

_{/ /} back to basics if (string (event .value) eq ' restore_pressed' ) then begin pixel_array [number_l - 1, *] = original_data [number_l - 1, *]

WIDGET_CONTROL, state.ymax, set_value = max (pixel_array ( (number_l - 1), *) ) y_max = max (pixel_array ( (number_l - 1), *)) WIDGET_CONTROL, state. ymin, set_value = min (pixel_array ( (number_l - 1), *) ) y_min = min (pixel_array ( (number_l - 1), *) )

,- _/ check for 10,-10 boundary... cheap... WIDGET_CONTRO , state.ymax, set_value = max([y_max, 10]) y_max = max([y_max, 10])

WIDGET_CONTROL, state. ymin, set_value = min([y_min, -10]) y_min = min([y_min, -10]) endif

,- ,- correct for zero-threshold anyways, it's cheap if ( (tl eq 0) OR (t2 eq 0)) then begin mess=WIDGET_MESSAGE(' Threshold cannot be 0!', /Error) WIDGET_CONTROL, state . threshold_one , SET_VALUE=thresholdl [number_l - l]

WIDGΞT_CONTROL, state . threshold_two, SET_VALUE=threshold2 [number_l - l] endif / ,- update the list of spikes time_array = where (yes_no_values (number_l - 1, *) eq 1) ,-this finds the location of the spikes for each cell number_spikes = n_elements (time_array)

,- ,- check for spike deletion if ( (time_array (0) ne -1) AND (spike_l ne 0)) then begin if ( (spike_l gt number_spikes) OR (spike_l It 1)) then begin mess=WIDGET_MESSAGE ( 'You entered an invalid spike number . ' )

WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0 endif else begin yes_no_values (nutnber_l - 1 , time_array ( spike_l - D ) = 0 draw_spikes, number_l _/that cell will now be displayed

WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0 endelse ,- ,- here is where all the plotting gets done endif else begin _/if no spikes are to be deleted draw_spikes, number_l

WIDGET_CONTROL, state . delete_spike, set_value = 0 spike_l = 0 endelse

,- ,- check for cell deletion if (string (event .value) eq 'delete_cell_pressed' ) then begin yes_no_values (number_l - 1,*) = 0 draw_spikes , number_l mess=WIDGET_MESSAGE (string ( ' Cell ' , byte (number_l) , ' has been successfully deleted. ' ) , /information) endif single_plot_defined=l /this is sent to

Multicell_event so the other programs now that single plots has been done end

Name : rms_spikes_wιdget

Description: a new, improved method of spike detection! Finds the mean and a number of standard deviations

; (generally 2... the Root Mean Square version of signal/noise handling) , and any points beyond this threshold are considered spikes .

History: finally created 3.26.99 by Rob.

4.20.99, Destroy widget upon completion. pro rms_spikes_widget common rms_spikes, rms_spikes_state rms_spikes_base=WIDGET_BASE (/COLUMN, title= ' Find Spikes- RMS') rms_spikes_button=WIDGET_BUTTON (rms_spikes_base , VALUE= ' Find Spikes', UVALUE=2) rms_spikes_type_bgroup=CW__BGROUP (rms_spikes_base , [ ' Positive Spikes' , 'Negative Spikes'], /ROW, /EXCLUSIVE, SET_VALUE=0) rms_spikes_threshold_field=CW_FIELD (rms_spikes_base, /RETURN_EVENTS , /FLOATING, TITLE= ' Threshold (Number of Std Devs) ',$ VALUΞ=2.0, UVALUE=0) rms_spikes_median_filter_window_size_field=CW_FIELD (rms_spikes_ base, /RETURN_EVENTS , /FLOATING, TITLE= ' Window Size for Median Filter' , $

VALUE=20.0, UVALUE=0) rms_spikes_text=WIDGET_TEXT (rms_spikes_base, VALUE='Set window size to 0 to prevent filtering...')

WIDGET_CONTROL, /realize, rms_spikes_base rms_spikes_state={ rms_spikes_base:rms_spikes_base, $ rms_spikes_button:rms_spikes_button, $ rms_spikes_type_bgroup : rms_spikes_type_bgroup, $ rms_spikes_threshold_field:rms_spikes_threshold_field, $ rms_spikes_median_filter_window_size_field : rms_spikes_median_fi lter_window_size_field}

WIDGET_CONTROL, WIDGET_INFO (rms_spikes_base, /CHILD), SET_UVALUE=rms_spikes_state xmanager, ' rms_spikes_widget ' , rms_spikes_base end

Name : rms_spιkes_widget_event Description:

History: Created February 25, 1999, by Rob.

Added the median-pass filter, 3.15.99, Rob.

3.26.99, Changed this to the event handler for the RMS spikes widget... a better, better world... _/ 3.26.99, Rewritten to be IDL 5.0 compatible... i.e. 'mean' and 'stddev' were replaced w/ general 'moment' _/ 4.20.99, Destroy widget upon completion. pro rms_spikes_widget_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, fi1e_name , $ total_frame_no, time_resolution, x_size, y_size, box_size common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined common with_choose_cells_com, thresholdl, threshold2, rms_threshold common rms_spikes, rms_spikes_state rms_state_stuff=WIDGET_INFO (event .handler, /Child) WIDGET_CONTROL, rms_state_stuff , GET_UVALUE=rms_spikes_state ,-find what is in /Child if (event. id eq rms_spikes_state.rms_spikes_button) then begin ,- ,- read interface data...

WIDGET_CONTROL, rms_spikes_state . rms_spikes_type_bgrou , GET_VALUE=pos_or_neg

WIDGET_CONTROL, rms_spikes_state . rms_spikes_threshold_field, GET_VALUE=threshold WIDGET_CONTROL , rms_spikes_state . rms_spikes_median_f ilter_window_size__f ield, GET_VALUE=window_s i ze for i = 0 , cell_no - 1 do begin rms_threshold [i] = threshold endfor if (window_size gt 0 ) then begin f ilter_median_subtractive , pixel_array, window_size , pixel_array end yes_no_values=intarr (cell no, frame_no)

,- ,- positive spikes is pos_or_neg = 0 if (pos_or_neg eq 0) then begin for cell_counter = 0, cell_no - 1 do begin cell_stats = moment (pixel_array [cell_counter, *] , sdev=cell_dff_stddev) cell_dff_mean = cell_stats [0] for frame_counter = 1, frame_no - 1 do begin if (pixel_array [cell_counter, frame_counter] gt (cell_dff_mean + (threshold * cell_dff_stddev) ) ) then begin yes_no_values [cell_counter, frame_counter] = 1 endif endfor endfor

,- ,- negative spikes is pos_or_neg = 1 endif else begin for cell_counter = 0, cell_no - 1 do begin cell_stats = moment (pixel_array [cell_counter, *] , sdev=cell_dff_stddev) cell_dff_mean = cell_stats [0] for frame_counter= 1, frame_no - 1 do begin if (pixel_array [cell_counter, frame_counter] It (cell_dff_mean - (threshold * cell_dff_stddev) ) ) then begin yes_no_values [cell_counter, frame_counter] = l endif endfor endfor endelse spikes_defined = 1 mess=WIDGET_MESSAGE( 'Spikes have been found!',

/INFORMATION)

WIDGET_CONTROL, rms_spikes_state . rms_spikes_base , /DESTROY endif end pro summed_spikes ^" common mother_com, pixel_array, yes_no_values , coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size base=WIDGET_BASE (/column, title= 'Overall slice behavior') sigma=CW_FIELD (base, title= ' Enter the value for sigma (smoothing factor):', VALUE=5, $

UVALUE=1 , /RETURN_EVENTS , /FLOATING) buttonl=widget_button (base , value= ' PLOT ' , uvalue= 'plot_pressed' )

WIDGET_CONTROL, base, /realize widget_control , base, set_uvalue=sigma xmanager, ' summed_spikes ' , base end

pro summed_spikes_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size

WIDGET_CONTROL, event. top, get_uvalue=sigma WIDGΞT_CONTROL, sigma, GΞT_VALUE=sigma window, 8 , title= ' Slice behavior¹

!except=0 summed_array=total (yes_no_values, 1) h=dblarr (frame_no) j =dblarr (frame_no) z=dblarr (frame_no) for t=0. , (frame_no -1), 1 do begin h[t]=0 j[t]=0 for s=0, (frame_no -1) do begin for i=l, Summed_array [s] do begin ,-the two below divisions are unnecessary h[t]=h[t] + (exp(l))^A(((-.5 * (s - t)^A2) ) / sigma^A2) / (Sqrt (2 * !pi * sigma^A2) ) endfor endfor for i=0, (frame_no-l) do begin j [t]=j [t] + (exp(l))^Λ(((-.5 * (i - t) 2) ) / sigma^A2) / (Sqrt (2 * !pi * sigma 2)) endfor z[t]=h[t]/j [t] endfor plot, f indgen (frame_no) *time_re solution, z, title= 'Plot of overall slice behavior', $ ^• xtitle = 'time in seconds', ytitle ='h(t) '

Unauthorized reproduction prohibited.

NAME:

TIFF_READ

PURPOSE :

Read 8-bit or 24-bit images in TIFF format.

CATEGORY :

Input/output .

CALLING SEQUENCE:

Result = TIFF_READ (Filename [,R, G, B] , ORDER = order, PLANARCONFIG = pc)

INPUTS :

Filename: A string containing the name of file to read. The default extension is ".TIF".

OUTPUTS :

TIFF_READ returns a byte array containing the image data. The dimensions of the result are the same as defined in the TIFF file: (Columns, Rows).

For TIFF images that are RGB interleaved by pixel, the output dimensions are (3, Cols, Rows).

For TIFF images that are RGB interleaved by image, on output Planarconfig is set to 2, and the result is the integer value zero. In this case, three separate images are returned in the R, G, and B output parameters.

OPTIONAL OUTPUTS:

R, G, B: Variables to hold the Red, Green, and Blue color vectors extracted from TIFF Class P, Palette Color images. For TIFF images that are RGB interleaved by image (Planarconfig returned as 2) the R, G, and B variables each hold an image with the dimensions (Columns, Rows) .

KEYWORDS :

The following keywords are used for output parameters only:

ORDER : The order parameter from the TIFF File. This parameter is returned as 0 for images written bottom to top, and 1 for images written top to bottom. If the Orientation parameter does not appear in the TIFF file, an order of 1 is returned .

_/ PLANARCONFIG: This parameter is returned as 1 for TIFF files that are Grayscale, Palette, or RGB color interleaved by pixel. This parameter is returned as 2 for RGB color TIFF files interleaved by image.

COMMON BLOCKS:

TIFF_COM. Only for internal use.

SIDE EFFECTS:

A file is read.

RESTRICTIONS :

Handles TIFF classes G, P, and R. One image per file.

EXAMPLE :

Read the file "my. tiff" in the current directory into the variable

; IMAGE, and save the color tables in the variables, R, G, and B by entering:

IMAGE = TIFF_READ("my.tiff", R, G, B)

To view the image, load the new color table and display the image by entering :

TVLCT, R, G, B TV, IMAGE

MODIFICATION HISTORY:

DMS, Written for VMS in 1985.

DMS, April, 1991. Rewrote and added class R and P images.

DMS, Jan, 1992. Fixed bug for images without a RowsPerStrip field.

DJC, Nov, 1993. Fixed doc header. DMS, Dec, 1994. Fixed bug with private tags. MWR, Mar, 1995. Fixed bug when opening non-existent file.

function tiff_long, a, i, len=len ; return longword(s) from array a(i) common tiff com, order, ifd, count on error, 2 _/Return to caller if an error occurs if n_elements (len) le 0 then len = 1 if len gt 1 then result = long (a, i, len) else result = long (a, i) if order then byteorder, result, /Iswap return, result end

function tiff_rational,a, i, len = len ,- return rational from array a (i) common tiff_com, order, ifd, count on_error,2 /Return to caller if an error occurs if n_elements (len) le 0 then len = 1 tmp = tiff_long(a, i, len = 2 * len) _/1st, cvt to longwords if len gt 1 then begin subs = lindgen(len) rslt = float (tmp (subs*2) ) / tmp (subs*2+l) endif else rslt = float (tmp (o) ) / tmp(l) return, rslt end function tiff_int,a, i, len=len /return unsigned long int from TIFF short int common tiff_com, order, ifd, count on_error,2 /Return to caller if an error occurs if n_elements (len) le 0 then len = 1 if len gt 1 then begin /Array? result = fix (a, i, len) if order then byteorder, result, /sswap result = long (result) if min (result) It 0 then begin /Convert to unsigned from signed 16bit negs = where (result It 0) result (negs) = 65535L + result (negs) endif endif else begin _/Scalar result = fix(a,i) if order then byteorder, result, /sswap if result It 0 then result = 65535L + result endelse return, result end function tiff_byte, a,i,len=len /return bytes from array a(i) common tiff_com, order, ifd, count on_error,2 /Return to caller if an error occurs if n_elements (len) le 0 then len = 1 if len gt 1 then result = a(i:i+len-l) else result = a(i) return, result end function tiff_read_field, index, tag, lun _/Return contents of field index

,- On output, tag = tiff tag index.

7 common tiff_com, order, ifd, count

on_error,2 /Return to caller if an error occurs TypeLen = [0, l, i, 2, 4, 8] /lengths of tiff types, 0 is null type for indexin ent = ifd (index * 12: index * 12 + 11) _/Extract the ifd tag = tiff_int (ent, 0) Tiff tag index typ = tiff_int (ent, 2) Tiff data type ent = tiff_long(ent, 4) # of elements nbytes = ent * TypeLen(typ) Size of tag field IF (nbytes GT 4) THEN BEGIN value size > 4 bytes ? offset = tiff_long(ent, 8) _/field has offset to value location

Point_Lun, lun, offset val = BytArr (nbytes) _/buffer will hold value (s)

Readu, lun, val

CASE typ OF _/Ignore bytes, as there is nothing to do

1 : i = 0 _/Dummy 2: val = String (val) ,-tiff ascii type

3: val = tiff_int (val, 0, len = ent) 4: val = tiff_long (val, 0, len = ent)

5: val = tiff_rational (val, 0, len = ent) ENDCASE ENDIF ELSE BEGIN _/Scalar...

CASE typ OF 1: val = ent (8)

2: val = string (ent (8 :8+ (cnt>l) -1) ) 3: val = tiff_int (ent, 8) 4: val = tiff_long (ent, 8)

ENDCASE ENDELSE return, val end

function tiff_read, file, r, g, b, order = ord, PlanarConfig = PC common tiff com, order, ifd, count

on_error,2 _/Return to caller if an error occurs openr, lun, file, error = i, /GET_LUN, /BLOCK if i It 0 then begin ,-OK? if keyword_set (lun) then free_lun,lun lun = -1 message, 'Unable to open file: ' + file endif hdr = bytarr ( 8 ) ,- Read the header readu, lun, hdr typ = string (hdr (0:1) ) _/Either MM or II if (typ ne 'MM') and (typ ne 'II') then begin message, 'TIFF_READ: File is not a Tiff file: ' + string(file) return, 0 endif order = typ eq 'MM' ,- 1 if Motorola 0 if Intel (LSB first or vax) endian = byte (1,0, 2) _/What endian is this? endian = endian (0) eq 0 ,-1 for big endian, 0 for little order = order xor endian ,-1 to swap...

_/ print, 'Tiff File: byte order=',typ, ', Version = ', tiff_int(hdr,2) offs = tiff_long(hdr, 4) _/Offset to IFD point_lun, lun, offs _/Read it a = bytarr(2) _/Entry count array readu, lun, a count = tiff_int (a, 0) _/count of entries

_/ print, count, ' directory entries' ifd = bytarr (count * 12) _/Array for IFD's readu, lun, ifd _/read it

,- Insert default values : compression = 1 bits_sample = 1 ord = 1 samples_pixel = IL pc = 1 photo = 1 rows_strip = 'fffffff 'xl _/Essentially infinity for i=0,count-l do begin _/Print each directory entry value = tiff_read_field (i, tag, lun) _/Get each parameter case tag of _/Decode the tag fields, other tags could be added

256 width = value

257 length = value

258 bits_sample = value

259 compression = value 262 Photo = value

273 StripOff = value

274 Ord = value

277 samples_pixel = long (value)

278 Rows_strip = value 279 Strip_bytes = value 284: PC = value 320: ColorMap = value else: value = 0 _/Throw it away endcase endfor

,- Do a cursory amount of checking: if bits_sample (0) ne 8 then $ message, 'TIFF_READ: bits_sample must be 8¹ if compression ne 1 then $ message, 'TIFF_READ: Images must be un-compressed' _/ if photo ge 4 then $

_/ message, 'TIFF_READ: Photometric Interpretation must be 0 to 3. ' if (pc eq 2) and (samples_pixel ne 3) then $ message, 'TIFF_READ: RGB data must have 3 SamplesPerPlane' nbytes = width * long (length) strips_image = (length + rows_strip -1) / rows_strip

if pc eq 1 then begin ,- lanar Config.... , simple if samples_pixel eq l then image = bytarr (width, length, /nozero) $ else image = bytarr (samples_pixel, width, length, /nozero) if strips_image eq 1 then begin _/Quick way? point_lun, lun, stripoff(O) ,-lst image data readu, lun, image _/Yes.... endif else begin ,-1 strip at a time.... tmp = Ob for i=0L, strips_image-l do begin point_lun, lun, stripoff (i) if n_elements (tmp) ne Strip_bytes (i) Then $ tmp = bytarr (Strip_bytes (i) , /nozero) readu, lun, tmp image (samples_pixel * width * i * rows_strip) = tmp endfor endelse if n_elements (ColorMap) gt 0 then begin _/Color map present? if n_elements (ColorMap) eq 768 then begin r = ishft (ColorMap (0 : 255) , -8) _/Remove and scale g = ishft (ColorMap (256: 511) , -8) b = ishft (ColorMap (512: 767) , -8) endif else message, ' IFF_READ : color map has wrong # of elements ' endif endif else begin ,-PC = 2, = interleaved by image

1 = 0 tmp = 0b for band = 0,2 do begin _/ Read each image image = bytarr (width, length, /nozero) for i=0L, strips_image-l do begin point_lun, lun, stripoff (1) if n_elements (tmp) ne Strip_bytes (1) then $ tmp = bytarr (strip_bytes (1) , /nozero) readu, lun, tmp image (width * i * rows_strip) = tmp

1 = 1 + 1 endfor _/Each strip case band of 0 : r = image 1 : g = image

2 : b = image endcase image = 0 endfor _/Each band endelse _/PC = 2 free_lun, lun return, image end

(c) , 1997 Copyright by Vikram Sheel Kumar Name: two_cells many times.pro ,- Synopsis: two_cells_many_times, no_iterations , least_no_of_matcb.es

/ Description: This procedure looks at the cells that fire together more than once. It then counts the number of times that pairs

/ of cells have fired more than once together. It also does this for random cells that have been created through the procedure

,- RANDOM_TEST. The random cells are tested for multiple hits through the procedure MULTIPLE_TEST_SIGNIFICANCΞ.

/ Explanation of Variables : no_iterations is the number of iterations the program will loop to create the random distribution.

/ least_no_of_match.es is the fewest number of times two cells have to spike together to begin counting the correlation.

_/ window_size is how far, to left and to right, of a spike we look for other spikes for two cells to be considered coactive

Created by VSK on June 1, 1997

Modified by VSK on June 2, 1997 to add least_no_of_matches

Modified by VSK on June 3, 1997 so now it doesn't look at the combinations of cells firing together more than x times but the number of times that

,- pairs of cells of the population pass the criteria and fire together more than or equal to lest_no_of_matches

Modified by VSK on June 18, to make it compatible with Mother4.0 Modified by VSK on July 1 to make it a part of Multicell Modified by VSK on July 7 to add get_lun and free_lun Modified by VSK on July 11, to add error checkers and replace get_lun and free_lun with open, close

,- Modified, Rob, 3.27.99, to add window size for connections

,- Modified, Rob, 4.7.99, to add 'mean distance between connections' statistic

_/ Modified, Rob, 4.14.99, speedup and cleanup.

,- Modified, Rob, 4.24.99, speedup- random test now ignores silent cells

,- Modified, 5.31.99, Rob, shifting to avoid wrapping pro two_cells_many_times, least_no_of_matches, no_iterations , window_size common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with_create_dist, cells_active ,-this contains the number of active cells forward_function multiple_test_significance common with two, filename connections_array = bytarr (cell_no, cell_no) yes_no_significance = yes_no_values

_/initializing the counters to zero true hits = 0

_/ creating the random distribution no_spikes = total (yes_no_values) random_distribution = intarr (no_iterations) _/this array holds one random number of matches for every iteration for t = 0, no_iterations - 1 do begin _/to create the random dist. we have to repeat the 2 steps no_iterations times make_random_data, seed, random_array, num_active_cells random_distribution(t) = count_random_hits_2_many (random_array, least_no_of_match.es , window_size, num_active_cells) endfor ; akes an array with a number for each iteration the number corresponds to the number of times 2 cells

^■ fire at least ' least_no_of_times ' in the random movie for the t-th iteration

_/printing some of the data filename=pickfile (/write, file= ( 'Two_Cells_' + strcompress (string (least_no_of_match.es) ) + '_Times ' ) ) if (filename eq ' ' ) then begin mess=WIDGET_MΞSSAGE( 'This data will not be saved!', /INFORMATION) filename= 'Two_Many.dat ' endif close, 1 openw, 1, filename printf, 1, 'Statistical Data:-' printf, 1, 'Number of times two cells must spike together to count as a hit : ' , least_no_of_matches printf, 1, 'Number of iterations: ', no_iterations printf, 1, 'Total number of active cells: ', cells_active _/ finding the actual number of matches distances_array = [0.0] for cell_l_counter = 0 , cell_no - 2 do begin for cell_2_counter = cell_l_counter + 1, cell_no - 1 do begin temp_hits = 0 temp_locations_cell_l = [-1] temp_locations_cell_2 = [-1] for window_counter = -window_size, window_size do begin temp_cell_l = intarr (frame_no + (2 * window_size) ) temp_cell_2 = intarr (frame_no + (2 * window_size) ) temp_cell_l (window_size : frame_no+window_size- 1) = yes_no_values (cell_l_counter, *) temp_cell_2 (window_size : frame_no+window_size- 1) = yes_no_values (cell_2_counter, *) temp_cell = (temp_cell_l * shift (temp_cell_2 , window_counter) ) _{/ / /} temp_cell = (yes_no_values (cell_l_counter,

*) * shift (yes_no_values (cell_2_counter, *) , window_counter) ) if ( (size (where (temp_cell) )) [0] eq 1) then begin temp_hits = (temp_hits + total (temp_cell) ) temp_locations_cell_l = [temp_locations_cell_l, where (temp_cell) - window_size] temp_locations_cell_2 = [temp_locations_cell_2 , where (temp_cell) - window_counter - window_size] endif endfor if (temp_hits ge least_no_of_matches) then begin

,- ,- rearrange location arrays- about the messiest code you can write, but it does the job:

,- ,- ...we want to first get rid of all -1 elements (meaning no paired spikes) , then sort the array, then (as sorting returns array subscripts

_{/ /} and not actual elements) we need to determine what the original values for the new subscripts are. If you're still confused,

/ / run this line of code in sections on test data, e.g. 'a = [-1,1,2,-1,4,8,-1,5] & a (sort (a (where (a ge 0) ) ) ) ' . real_locations_cell_l = (temp_locations_cell_l (where (temp_locations_cell_l ge 0) ) ) (sort (temp_locations_cell_l (where (temp_locations_cell_l ge 0) ) ) ) real_locations_cell_2 = (temp_locations_cell_2 (where (temp_locations_cell_2 ge 0))) (sort (temp_locations_cell_2 (where (temp_locations_cell_2 ge 0))))

,- ,- compute connection distance connection_distance = sqrt ( ( (double (location (cell_l_counter) . coord [0] ) - double (location (cell_2_counter) .coord [0] ) ) 2) + ( (double (location (cell_l_counter) . coord [1] ) - double (location (cell_2_counter) .coord [1] ) ) ^A 2)) distances_array = [distances_array, connection_distance] connections_array (cell_l_counter, cell_2_counter) = 1 true_hits = true_hits + 1 ; ,- dump to file printf, 1, 'Cell number one and location: strcompress (cell_l_counter + 1) , location (cell_l_counter) .coord printf, 1, 'Frames cell one fires in:' printf, 1, real_locations_cell_l printf, 1, 'Cell number two and location: strcompress (cell_2_counter + 1) , location (cell_2_counter) .coord printf, 1, 'Frames cell two fires in:' printf, 1, real_locations_cell_2 printf, 1, ' ' / / update 'significance array' which is used to draw the raster plot for frame_counter = 0 , temp_hits - 1 do begin yes_no_signif icance (cell_l_counter, real_locations_cell_ι [frame_counter] ) =TΞMPORARY (yes_no_signif icance (c ell_l_counter, real_locations_cell_l [f rame_counter] ) ) + 1 yes_no_significance (cell_2_counter, real_locations_cell_2 [frame_counter] ) =TEMPORARY(yes_no_significance (c ell_2_counter, real_locations_cell_2 [frame_counter] ) ) + 1 endfor endif endfor endfor

_/draw raster plot before it checks to see if the random distribution has a variance as the raster is independant of statistics draw_significance_raster, yes_no_significance

_/check to see if event occurs randomly and thus if a p value can be calculated. xmax=max (random_distribution) xmin=min (random_distribution) if (xmax eq xmin) then begin mess=WIDGET_MESSAGE( 'Random distribution has variance of zero. Try again with greater number of iterations', /error) printf, 1, 'ERROR!' printf, 1, 'Moment undefined for random distribution with variance zero' free_lun, 1 close, 1 p_value=-l endif else begin

_/making the histogram for the random distribution !p.multi=0 bin_size=l hist=histogram(random_distribution, binsize=l, min=0, max= (2*xmax) +1) _/plot a histogram of the distribution window, /free, title='Two Many Distribution' plot, hist, xtitle= 'number of hits', ytitle=' frequency ' y2=total (random_distribution eq true_hits) _/for drawing the blue line for the real data, we have to find the height of that line in the histogram if y2 le 0 then begin _/if there are no random values equal to the true data, draw a line of height one y2=l endif plots, [true_hits, true_hits] , [0, y2] , color=12 ,-this draws a line where the actual number of matches lies

_/calculating the number of spikes per cell per second spikes_per_cell_per_second= (no_spikes/cells_active) / (frame_no*t ime_resolution)

_/draw the correlation map connections_array= (connections_array) +TRANSPOSΞ (connections_arr ay) _/this is to make the connections_array symmetric. Because if cell 1 is connected to 2 , 2 is also connected to 1. correl_map_plane needs the symmetric array window, /free, xsize=x_size, ysize=y_size, title= ' Correl Map for Two Many¹ correl_map_image_plane, connections_array _/find the p value, standard deviation etc. stats=moment (random_distribution, sdev=sdev) no_points_right=total (random_distribution ge true_hits)_. p_value=no_points_right/no_iterations _/print the data to a file printf, 1, 'Total number of frames: ', strcompress (frame_no) printf, 1, 'Total number of spikes: ', strcompress (no_spikes) printf, 1, 'Mean expected matches: ', strcompress (stats (0) ) printf, 1, 'Variance: ', strcompress (stats (1) ) printf, 1, 'Standard deviation: ', strcompress (sdev) printf, 1, 'Actual matches: ', strcompress (true_hits) printf, 1, 'Actual/expected: ', strcompress (true_hits/stats (0) ) printf, 1, 'Standard error for ratio: ', strcompress (sdev/stats (0) ) printf, 1, 'Spike firing rate: ', strcompress (spikes_per_cell_per_second) if (n_elements (distances_array) gt 2) then begin distance_stats = moment (distances_array [1 :n_elements (distances_array) -1] , sdev=sdev_distances_array) printf, 1, 'Mean connection distance: ', strcompress (distance_stats [0] ) printf, 1, 'Standard deviation: ', strcompress (sdev_distances_array) endif else begin if (n_elements (distances_array) gt 1) then begin printf, 1, 'Connection distance: ', strcompress (distances_array [1] ) endif endelse printf, 1, 'Significance p-value: ', strcompress (p_value) close, 1 endelse end

^~; (c) , 1997 Copyright by Vikram Sheel Kumar _/ Name : two_many_widget

/ Synopsis: two_many_widget,

/ Description: Using this widget the user can test the significance of two cells firing many times in

; the movie to analyze. The user will have to input the specifications such as the number of times

/ that two cells fire together to be taken as the min. criterion for a hit. / Modification History: modified 3.27.99, Rob, to add window size to cell-connection finding

/ modified, 4.21.99, Rob, to only do the test when the button is pressed and to init fields to default values pro two_many_widget common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, fi1e_name , $ total_frame_no , time_resolution, x_size, y_size, box_size two_many_base=WIDGET_BASE (/COLUMN, title= ' Significance of two cells firing together many times') fieldl=CW_FIELD(two_many_base, /RETURN_E ENTS, /INTEGER, TITLE='Two cells should fire together at least this many times:', VALUE=2 , UVALUE=0 ) field2=CW_FIELD(two_many_base, /RETURN_EVΞNTS , /INTEGER, TITLE= 'Number of iterations:', VALUE=1000, UVALUE=0) field3=CW_FIELD(two_many_base, /RΞTURN_EVENTS , /INTEGER, TITLE= ' Window size for hits:', VALUE=0, UVALUE=0) buttonl=WIDGET_BUTTON(two_many_base, VALUE= ' Find significance', UVALUE=2)

WIDGET_CONTROL, /realize, two_many_base two_many_state= {fieldl : fieldl, field2 : field2 , field3 :field3 , buttonl : buttonl}

WIDGET_CONTROL, WIDGET_INFO (two_many_base, /CHILD), SET_UVALUE=two__many_state xmanager, ' two_many_widget ' , two_many_base end

©Copyright Synapse Innovation Inc . Written by Vikram Sheel Kumar NAME : two_many_widget_event

SYNOPSIS: two_many_widget_event , event (called by xmanager)

DESCRIPTION: This is the event handler for the two_many_widget .

HISTORY: Created on June 18, 1997 by Vikram Sheel Kumar

Modified on July 26, 1997 by VSK to check for size of file before xdisplaying

Modified, 3.27.99, Rob, to add window size for connection finding pro two_many_widget_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common with_two, filename two_many_state_info=WIDGΞT_INFO (event . handler, /CHILD) WIDGET_CONTROL, two_many_state_info, GET_UVALUE=two_many_state WIDGΞT_CONTROL, two_many_state . fieldl , GET_VALUE=least_no_of_match.es WIDGET_CONTROL, two_many_state . field2 ,

GET_VALUE=no_of_iterations

WIDGET_CONTROL, two_many_state . field3 , GET_VALUE=window_size if (event. id eq two_many_state .buttonl) then begin le 1) !',

endif else begin two_cells_many_times , least_no_of_matches , no_of_iterations, window__size WIDGET_CONTROL, event. top, /hourglass close, 1 openr, 1, filename stat_stuff=fstat (1) file_size=stat_stuff . size close, 1 if (file_size gt 8112) then begin mess=WIDGET_MESSAGE ( 'File is too large to display through a widget. Open it manually. If you did not save it, it is named Two_Many.dat', /INFORMATION) endif else begin xdisplayfile, filename, title = "Statistical Data for Two Many", group = event.top, width = 75, height = 50 endelse endelse endif end

,- Widget for creating correlation coefficients matrix and drawing general correlation maps .

7 ,- Modified, 4.19.99, Rob, to change a couple of names... pro widget_analyze common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, $ total frame_no, time resolution, x_size, y size, box_size

basel=WIDGET_BASE (/COLUMN, title= ' General Slice Correlation') type={cw_pdmenu_s, flags:0, name:'¹} nitty_gritty =[ { cw_pdmenu_s, 0, 'Calculate Correlation Coeff Matrix' }, $

{ cw_pdmenu_s, 1, 'Build Correlation Map'}, $ { cw_pdmenu_s , 0 , 'On Image ' } , $ { cw_pdmenu_s, 2, 'Spatial'}] pull_down2=CW_PDMΞNU(basel, nitty_gritty, /RETURN FULL NAME)

WIDGET_CONTROL, basel, /realize xmanager, 'widget_analyze ' , basel end

(c) Copyright, 1996 by Neuron Stimulation Inc. Name : widget_analyze_event Synopsis: widget_analyze_event, event

Description.- This procedure is called when you press the 'Analyze' button on the main menu. / It allows you to calculate the correlation coefficient (WHICH MUST

BE THE FIRST STEP) and then

_/ draw a correlation map on an image or spatial .

7

,- Modified on Nov 19, 1996

,- Modified, 4.19.99, Rob, to change a couple of names...

EVENT MANAGER FOR ANALYSIS

WIDGET pro widget_analyze_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name,$ total_frame_no, time_resolution, x_size, y_size, box_size common flags, cells_defined, spikes_defined, single_plot_defined, correl_coef_defined

CASE (event.value) of

'Calculate Correlation Coeff Matrix' : BEGIN choose_correl END

'Build Correlation Map.On Image': BEGIN if (correl_coef_defined ne 1) then begin mess=WIDGET_MESSAGE ( 'You have to first calculate the correlation coefficients ! ' , /Error) endif else begin correl_map_image endelse END

'Build Correlation Map. Spatial¹ : BEGIN if (correl_coef_defined ne 1) then begin mess=WIDGET_MΞSSAGΞ( 'You have to first calculate the correlation coefficients ! ' , /Error) endif else begin base5=WIDGET_BASΞ (/COLUMN, title= ' Spatial Correlation Map ' ) draw7=WIDGET_DRAW(base5, xsize=x_size, ysize=y_size)

WIDGET_CONTROL, /realize, base5 correl_map_image_plane, coef endelse END

ENDCASE end APPENDIX B

_/ to find argument of the maximum value of the array ,- used by hidden_markov_event .pro function arg_max, x,y,j,t,D,A

Z=dblarr (y+1) for i=x,y do begin if ( A[i,j] eq 0) then begin z [i]=-10D307 endif else begin z[i] = D[t-l,i]+alog(A[i, j]) endelse endfor m=max (Z, k) return, k end

; function used by HMM to find the probability of observation vector of

_/ length cell_no, at time t, given the array of Poisson lambdas B, and state_no function B_prob, state_no, observation, t, B, cell no

_/ ! except=2 s=double(l) for i=0, cell_no -1 do begin observation (i, t) =double (observation (i,t) ) if (B [state_no, i] eq 0) then begin if (observation (i, t) eq 0) then begin s=s * 1. endif else begin s=s * 0 endelse endif else begin s=s* ( (exp (-B [state_no,i] ) ) * ( (B [state no, i] ) ^λ (observation (i, t) ) )

/ factorial (observation (i, t) ) )

,-s=s - B [state_no, i] + (observation (i, t) ) *alog (B [state_no, i] ) - alog (factorial (observation (i, t) ) ) endelse endfor

_/ s=exp (s) return, s end pro choose__init_par

common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common markov, n, m, A, B, P, observation, f , g, d, Fi , Q, back, xsi , Gamma, n_i , a_i , b_i , p_i , c , _s

, t_max common indexes , i_a, i_b, i_p i_a=0 i_b=0 i_p=0 Values= [ 'random' , 'uniform' ]

Values2= [ ' random' , 'uniform' , ' averaged over intervals ' ] base=widget_base (/column, title=' Initial parameters generation')

Ainit=widget_droplist (base, value=values, uvalue=0, title= ' Choose initial state transition prob. A: ')

Binit=widget_droplist (base, value=values2, uvalue=l, title= ' Choose initial state characteristcs B: ')

Pinit=widget_droplist (base, value=values, uvalue=2, title= ' Choose initial state prob. P: ') Apply=widget_button(base, value= 'Apply ' , uvalue=3) widget_control , base, /realize xmanager, ' choose_init_par ' , base end pro choose_init_par_event, event

common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common markov, n,m, A,B, P, observation, f , g,d, Fi, Q, back, xsi, Gamma, n_i, a_i,b_i,p_i, c, w_s

, t_max common indexes, i_a,i_b,i_p

WIDGET_CONTROL, event . id, GET_UVALUE=uval

IF (TAG_NAMES (event, /STRUCTURE_NAME) EQ ' WIDGET_DROPLIST' ) $ THEN BEGIN Case uval of 0 : begin CASE event ^'. index OF

0 : begin

_/print, 'A random' i_a=i_a+l _/Random distribution for i=0, (n-1) do begin for j=0, (n-1) do begin

A[i, j] =randomu (seed) endfor endfor zl=dblarr (n) zl=total (A, 2) for i=0, (n-1) do begin for j=0, (n-1) do begin A[i, j]=A[i, j]/zl[i] _/????? endfor endfor end

1 : begin

_/print, 'A uniform' i_a=i_a+l

_/Uniform distribution for i=0, (n-1) do begin for j=0, (n-1) do begin

A[i,j]=l/n endfor endfor end

ENDCASE end : begin

CASE event . index OF

0 : begin _/print, 'B random' i_b=i_b+l for i=0, (n-1) do begin for j=0, (m-1) do begin

B [i, j ] =Randomu(seed) _/Random endfor endfor zl=total (b,2) for i=0, (n-1) do begin for j=0, (m-1) do begin B[i, j]=B[i, j]/zl[i] _/think endfor endfor end 1: begin

_/print, 'B uniform¹ i_b=i_b+l

_/ uniform for i=0, (n-1) do begin for j=0, (m-1) do begin

,-mm=m .

B[i, j]=l/float(m) endfor endfor end

2 : begin _/ average over uniform segments i_b=i_b+l for j=0, m-1 do begin for i=0, n-1 do begin s=0 s i=0 for k=f ix(i* ( (frame_no - w_s + 1) /n) ) , f ix( (i+l) * ( (frame_no -w_s + l)/n)-l) do begin _/change frame_no s=s+observation [ j , ] s_i=s_i+l ,- print , k endfor

_/print, 'dupa' s=double (s) s_i= double (s_i) B[i, j]=double( s / s_i) endfor endfor end ENDCASE end begin CASE event . index OF

0 : begin

_/print , ' P random ' i_p=i_P+l for i=0, (n-1) do begin P (i) =randomu(seed) _/random endfor z2=total (p) for i=0, (n-1) do begin

P(i)=P(i)/z2 endfor end

1 : begin

_/print, 'P uniform' i_p=i_p+l

_/ uniform

P(0)=1 for i=l, (n-1) do begin

P(i)=0 endfor end

ENDCASE end endcase ENDIF

if (uval eq 3) then begin print, 'Apply pressed' /print, ' i_a= ' , i_a, ' , i_b=',i_b,', i_p= ' , i_p if (i_a eq 0) then begin _/make random _/print, 'A random' for i=0, (n-1) do begin for j=0, (n-1) do begin

A[i, j] =randomu (seed) endfor endfor zl=dblarr (n) zl=total(A,2) for i=0, (n-1) do begin for j=0, (n-1) do begin

A[i, j]=A[i, j]/zl[i] ,-????? endfor endfor endif if (i_b eq 0) then begin _/ make random

_/print, 'B random¹ for i=0, (n-1) do begin for j=0, (m-1) do begin

B [i, j] =Randomu (seed) endfor endfor zl=total (b,2) for i=0, (n-1) do begin for j=0, (m-l) do begin B[i, j]=B[i, j]/zl[i] _/think endfor endfor endif if (i_p eq 0) then begin _/ make random _/print, 'P random' for i=0, (n-1) do begin

P (i) =randomu(seed) endfor z2=total (p) for i=0, (n-1) do begin

P(i)=P(i)/z2 endfor endif print, 'These are initial state transition probabilities, A' print, A print, 'space' print, 'These are initial probabilities of observing symbol m at state n, B' print, b print, 'space' print, 'This is the initial state distribution, P' print, p print, 'space' n_i=n _/for save function

A_i=A

B_i=B

P_i=P widget_control , event. top, /destroy endif end pro hidden_markov common mother_com, pixel_array, yes_no_values , coef , location, cell_no , frame_no , file_name , direct , $ total_frame_no , time_resolution, x_size , y_size , box_size common markov, n , m, A, B , P , observation, f , g , d , Fi , Q , back, xsi , Gamma , n_i , a_i , b_i , p_i , c , w_s

, t_max common markov2 , yes_no_values3 , yes_no_temp , cell_no_temp , frame_no_temp

this is the definition of a test input with 400 frames and 4 states frame_no=400 observation=intarr (frame_no) _/observation sequence for i=0,99 do begin x=randomn (seed,binomial= [1, .3]) if (x eq 1) then begin observation [i] =1 endif else begin observation [i] =2 endelse endfor for 1=100,199 do begin x=randomn(seed,binomial= [1, .6] ) if (x eq 1) then begin observatio [i] =1 endif else begin observation [i] =3 endelse endfor for i=200,299 do begin x=randomn(seed,binomial= [1, .1]) if (x eq 1) then begin observation [i] =2 endif else begin observation [i] =0 endelse endfor for i=300,399 do begin x=randomn(seed,binomial= [1, .5] ) if (x eq 1) then begin y=randomn (seed,binomial= [1, .4] ) if (y eq 1) then begin observation [i] =1 endif else begin observation [i] =2 endelse endif else begin y=randomn(seed,binomial= [1, .4] ) if (y eq 1) then begin observation [i] =0 endif else begin observation [i] =3 endelse endelse endfor m=double (4) /this transforms yes_no_values to observation seq., information is lost because only one cell per time unit is allowed to spike observation=intarr (frame_no) /observation sequence for i=0, (frame_no -1) do begin observation [i] =0 for j=0, (cell_no -1) do begin if (yes_no_values [j , i] eq 1) then begin observation [i] =j+l /needs change endif endfor endfor m=double (cell_no + 1) / number of observation symbols (+1 if neither cell fires - 0)

/this transforms yes_no_values to observation by adding everything - creates artificially many zeros

,-m=double (cell_no + 1) / number of observation symbols (+1 if neither cell fires - 0) observationl=intarr (frame_no * cell_no) /observation sequence z=0 for j=0, (frame_no -1) do begin for i=0, (cell_no -1) do begin if (yes_no_values [i, j] eq 1) then begin observationl [z] =i+l endif else begin observationl [z] =0 endelse z=z+l endfor endfor frame no=frame no*cell no

,-this creates observation from yes-no_values by adding everything but only if cells fire simulatneously

,-m=double (cell_no + 1) / number of observation symbols (+1 if neither cell fires - 0) observation_2=intarr (frame_no * cell_no) /observation sequence z=0 for j=0, (frame_no -1) do begin z2=0 /for i=0, (cell_no -1) do begin if (yes no_values [i, j] eq 1) then begin observation_2 [z] =i+l z=z+l _. z2=l ,- endif else begin observation_2 [z] =0 endelse endfor if (z2 eq 0) then begin z=z+l endif endfor

,-observation2=intarr (z)

,-for i=0, (z-1) do begin

/ observation [i] =observation_2 [i]

,- endfor

_/ frame no=z

_/ observation=observation2

test=fltarr(2,3) test [0,0] =1.1 test [0,1] =1.2 test [0,2]=1.3 test [1,0] =2.1 test [1,1] =2.2 test [1,2] =2.3 suml=fltarr (5) suml=total (test, 1) sum2=fltarr (5) sum2=total (test, 2) print, test print, 'this is sum print, su l print, 'this is sum print, sum2

base=widget_base (/column, title= 'Hidden Markov Modelling') stuff= {cw_pdmenu_s , flags :0, name:¹'} details=[{ cw_pdmenu_s, 1, 'Create input data' }, $

{ cw_pdmenu_ -S, 0 ' from observation sequence ' } , $

{ cw_pdmenu _^s' 0 ' from shuffled in time observation sequence ' } . $

{ cw_pdmenu ³* 0 'from Monte Carlo simulated observation seq

(space) ' },$

{ cw_pdmenu _^Ξ' 2 'test' },$

{ cw_pdmenu_ _^s' 0 ' Record y_n_values ' } , $ { cw_pdmenu_s, 0, 'Retrieve y_n_values ' } ] pull_down=cw_pdmenu (base, details, /return_name, UVALUE=4) buttonl=cw_bgroup(base, /row, ['Choose initial parameters.', 'Save init. pars.', 'Load init. pars.', 'Save final pars.', 'Load final pars.'], /return name, UVALUE=0) state_no=CW_FIELD (base, title= ' Enter the number of states you would like to use in the model:', VALUE=4, $ UVALUE=5, /RETURN_EVΞNTS, /Floating) button2=widget_button(base, value= 'Estimate model parameters.¹, uvalue=3 ) iteration_no=CW_FIELD (base , title= ' Enter the number of iterations you would like to perform: ' , VALUΞ=20 , $

UVALUE=4, /RΞTURN_ΞVENTS , /FLOATING) button3=widget_button(base, value='Find the most probable hidden state sequence.', uvalue=2) button4=widget_button(base, value=' Compute P (observation | model parameters) ' , uvalue=l) button5=widget_button(base, value='Find cross correlations within states', uvalue=6) correlation_state=CW_FIELD (base, title= ' In which state do you want to find cross-correlations ?:', VALUE=0, $

UVALUE=7, /RETURN_EVENTS, /FLOATING) state= { state_no : state_no, iteration_no : iteration_no, correlation_state : correlation_state} widget_control , base, set_uvalue=state

WIDGET_CONTROL, state_no, GET_VALUE=n widget_control , base, /realize xmanager, ' hidden_markov ' , base end

pro hidden_markov_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no , time_resolution, x_size, y_size, box_size common markov, n,m, A,B, P, observation, f ,g,d, Fi,Q, back, xsi, Gamma, n_i, a_i,b_i,p_i, c, w_s , t_max common markov2 , yes_no_values3 , yes_no_temp , cell_no_temp , frame_no_temp

_/ ! except=2 n_old=n

WIDGET_CONTROL, event. id, GET_UVALUE=uval

Widget_control , event. top, get_uvalue=state

WIDGET_CONTROL, state . iteration_no, GET_VALUE=iteration_no

WIDGET_CONTROL, state . state_no, GET_VALUE=n WIDGET CONTROL, state . correlation state, GET VALUE=correlation state

CASE uval OF 4 : Begin

Case event.value OF

'from observation sequence': BEGIN print, 'these are y_n_values ' print, yes_no_values print, ' ' m=cell_no w_s= 30. / size of sliding window, step is equal to 1 observation=intarr (cell_no, FRAME_NO - w_s +1) for i=0, (frame_no - w_s ) do begin for k=0, (cell_no -1) do begin s=0 if (i It (frame_no - w_s +1) ) then begin for j=i, (i+w_s -1) do begin s=s+yes_no_values [k, j ] endfor observation [k, i] =s endif endfor endfor print, 'Input from observation sequence created' end

'from shuffled in time observation sequence': BEGIN

yes_no_values2=intarr (cell_no, frame_no) shuffle=intarr (cell_no) for i=0, cell_no -1 do begin shuffie [i] = (frame_no / 2) * randomn (seed, uniform=l) endfor for i=0, cell_no -1 do begin for j=0, frame_no -1 - shuffie [i] do begin yes_no_values2 [i , j ] =yes_no_values [i , j +shuffle [i] ] endfor for j=(frame_no - shuffie [i] ) , (frame_no -1) do begin yes_no_values2 [i,j ] =yes_no_values [i, shuffie [i] - frame_no + j] endfor endfor

/print, 'These are yes and no values:' /print, yes_no_values2 /print, ' ' ,-yes_no_values=yes_no_values2 m=cell_no w_s= 30. / size of sliding window, step is equal to 1 observation=intarr (cell_no, FRAME_NO - w_s +1) for i=0, (frame_no - w_s ) do begin for k=0, (cell_no -1) do begin s=0 if(i It (frame_no - w_s +1)) then begin for j=i, (i+w_s -1) do begin s=s+yes_no_values [k, j ] endfor observation [k, i] =s endif endfor endfor print, 'Input from shuffled in time observation sequence created' end 'from Monte Carlo simulated observation seq (space) ': BEGIN _/this is copied from vikram' s make_random_data.pro

_/find the number of spikes each cell of the true data has no_spikes=intarr (cell_no) ,-this is the array with the number of spikes for each cell for cell=0, cell_no-l do begin no_spikes (cell) =total (yes_no_values (cell, *) ) endfor cells_active=total (no_spikes gt 0)

_/create a random binary spike train. First, numbers for locations of spikes are generated randomnly making

,- sure that no number is repeated more than once. Next, a binary spike train is created using these locations

,-for the position of spikes. Each cell has the same number of spikes as the true data had. random_array=intarr (cell_no, frame_no) for cell=0, cell_no-l do begin index=no_spikes (cell) -1 _/using a random number generator for a uniform distribution

_/here we see if there are any spikes in the original cell of the data if index ge 0 then begin repeat begin temp=randomu(my_seed, no_spikes (cell) ) _/random numbers temp=fi (frame_no*temp)

_/this has the location of the unique elements in the random row of data

_/this is the test to make sure no numbers are repeated in the data unique= (uniq (temp (sort (temp) ) ) ) _/# of unique elements numberl=n_elements (unique)

_/total number of elements number2=n_elements (temp) endrep until numberl eq number2 random_array (cell, temp-l)=l

_/put spikes at the random locations endif endfor yes_no_values=random_array _/print, 'this are y_n values' _/print, yes_no_values m=cell_no w_s= 30. _/ size of sliding window, step is equal to 1 observation=intarr (cell_no, FRAME_NO - w_s +1) for i=0, (frame_no - w_s ) do begin for k=0, (cell_no -1) do begin s=o if (i It (frame_no - w_s +1) ) then begin for j =i , (i+w_s -1) do begin s=s+yes_no_values [k, j ] endfor observation [k, i] =s endif endfor endfor print, 'Input from monte carlo simulated data created' end

' test ' : BEGIN

,- this is going to be a test print, 'Frames 1-100 should be in state 0 with paramters : [ 3,2,8,1] ' print, 'Frames 101-200 should be in state 1 with parameters: [ 3,4,8,4] ' print, 'Frames 201-300 should be in state 2 with parameters: [ 3,3,6,0] ' print, 'Frames 301-400 should be in state 3 with parameters: [11,3,6,4] •

m=4. cell_no=4. w_s=ll frame_no=410 observationsintarr (cell_no, FRAME_NO - w_s +1) for i=0, (99 ) do begin observation 0, i] =randomn (seed, poisson=3) observation 1, i] =randomn(seed, poisson=2) observation 2, i] =randomn(seed, poisson=8) observation 3, i] =randomn(seed, poisson=l) endfor for i=100, 199 ) do begin observation 0, i] =randomn(seed, poisson=3) observation 1, i] =randomn(seed, poisson=4) observation 2, i] =randomn(seed, poisson=8) observation 3, i] =randomn(seed, poisson=4) endfor for i=200, 299 ) do begin observation 0, i] =randomn(seed, poisson=3) observation 1, i] =randomn(seed, poisson=3) observation 2, i] =randomn(seed, poisson=6) observation 3,i]=0 endfor for i=300, 399 ) do begin observation 0, i] =randomn(seed, poisson=ll) observation 1, i] =randomn(seed, poisson=3) observation 2, i] =randomn(seed, poisson=6) observation 3 , i] =randomn(seed, poisson=4) endfor print , ' test input created ' end 'Record y_n_values ' : begin yes_no_temp=yes_no_values cell_no_temp=cell_no frame_no_temp=frame_no print, 'Yes and no values recorded' end

'Retrieve y_n_values ' : begin cell_no=cell_no_temp frame_no=frame_no_temp yes_no_values=intarr (cell_no, frame_no) yes_no_values=yes_no_temp print, 'Yes and no values retrieved' end

' test2 ' :begin for i=0,4 do begin print, yes_no_values [i, 7] endfor end else: endcase _/print, 'This is the input:' ,-print, observation _/print, ' ' t_max=frame_no - w_s +1 n=double (n)

A=dblarr (n,n) B=dblarr (n,m) P=dblarr (n) end

0 : Begin

Case event.value OF

' Choose initial parameters . ' : BEGIN

if (n_old ne n) then begin n=double (n)

A=dblarr (n,n)

B=dblarr (n,m)

P=dblarr (n) endif choose_init_par end

' Save init . pars . ' : BEGIN data_file=pickfile (/write, file='HM_IVars_exp' , title= 'Create the saved variables file') if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGE ( 'No saved variables file specified', /INFORMATION) endif else begin

Save, FileName=data_file, A_i, B_i, P_i, n_i mess=WIDGET_MESSAGE( ' Initial data from this experiment saved' , /INFORMATION) endelse end

'Load init. pars.': BEGIN data_file=pickfile (/read, title= 'Select the saved variables file¹, GΞT_PATH=filepathl) if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGE ( 'No saved variables file specified', /INFORMATION) endif else begin restore, data_file mess=WIDGET_MESSAGE ( 'Previous initial A, B, P, and n variables loaded', /INFORMATION) endelse data_file=0 a=a_i b=b_i p=p_i n=n_i print, 'This is initial A' print, A print, ' ' print, 'This is initial B' print, B print, ' ' print, 'this is initial p' print, p print, ' '

WIDGET CONTROL, state. state no, SET VALUE=n end

' Save final pars . ' : BEGIN data_file=pickfile (/write, file='HM_FVars_exp' , title= ' Create the saved variables file ' ) if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGE( 'No saved variables file specified', /INFORMATION) endif else begin

Save, FileName=data_file, A, B, P, n mess=WIDGET_MESSAGE( 'Data from this experiment saved', /INFORMATION) endelse end

' Load final pars . ' .- BEGIN data_file=pickfile (/read, title= ' Select the saved variables file', GET_PATH=filepathl) if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGE( 'No saved variables file specified', /INFORMATION) endif else begin restore, data_file mess=WIDGET_MESSAGE( 'A, B, P, and n variables loaded', /INFORMATION) endelse data_file=0 print, 'This is A' print, A print, ' ' print, 'This is B' print, B print, ' ' print, 'this is p' print, p print, ' '

WIDGET_CONTROL, state . state_no, SET_VALUE=n end

else : endcase end

1 : Begin

/The Log probability of the observation given the model g=double (0) for i=0, n-1 do begin g=g+f (t_max-l, i) _/ endfor for t=0,t_max -1 do begin g=g - alog(c(t)) endfor print, 'This is the Log probability of the observation sequence given the model ' print, g print, ' ' end

,-2 :begin _/ Viterbi with no scaling _/print, 'button 2 pressed'

,- This is an implementation of the viterbi algorithm, the pupose of _/ viterbi algorithm is to find the single most probable state sequence given / the observation sequence and model parameters.

_/ !except=2

,-D=dblarr (t_max,n) _/highest probability along a single state path of the observation sequence ,-Fi=dblarr (t_max,n) / the backtracking step to retrieve the hidden states initiation for i=0, n-1 do begin

D [0, i] =p [i] *B_prob (i, observation, o)

Fi [0,i]=0 endfor recursion

Z=dblarr (n) for t=l, t_max-l do begin for j=0, n-1 do begin for i=0, n-1 do begin

Z[i]=D[t-l,i]*A[i,j] endfor

D [t, j] = (max(Z) ) *B [j , (observation (t) ) ]

Fi [t, j] =arg_max(0,n-l, j , t,D, A) ,- i changes from 0 to n-1 endfor endfor termination for i=0, n-1 do begin Z [i] =D [t_max-l,i] endfor

Prob_path=max(Z,k) ,- probability of the given path Q=intarr (t_max) _/ Q is the hidden state sequence Q [t_max -l]=k for t=t_max-2, 0, -1 do begin Q[t]=Fi[t+l, (Q[t+1])] endfor print, 'This is D' print, D print, 'this is Fi ' print, Fi print, 'This is Q' print, Q print, ' ' print, 'This is the corresponding probability of the hidden state sequence given the observation data: ' print, Prob_path print, ' ' end 2 .-begin

,-print, 'button 2 pressed'

_/ This is an implementation of the viterbi algorithm, the pupose of ; viterbi algorithm is to find the single most probable state sequence given

_/ the observation sequence and model parameters . here i use logarithmic scaling

D=dblarr (t_max,n) /highest probability along a single state path of the observation sequence

Fi=dblarr (t_max,n) _/ the backtracking step to retrieve the hidden states

_/ initiation for i=0, n-1 do begin if( (P(i) eq 0) or (B_prob (i, observation, 0,B, cell_no) eq 0) ) then begin

D[0,i]=-10D307

Fi [0,i]=0 endif else begin

D [0, i] =alog(P (i) ) +alog(B_prob (i, observation, 0,B, cell_no) ) Fi[0,i]=0 endelse endfor

_/recursion Z=dblarr(n) for t=l, t_max-l do begin for j=0, n-1 do begin for i=0, n-1 do begin if ( A [i , j ] eq 0) then begin Z [i] =-10D307 endif else begin

Z [i] =D [t-l , i] +alog (A [i , j ] ) endelse endfor if ( B_prob (j , observation, t,B, cell_no) eq 0) then begin

D[t, j]=-10D307

Fi [t, j] =arg_max(0,n-l, j , t,D,A) ,- i changes from 0 to n-1 endif else begin

D [t, j] = (max(Z) ) +alog(B_prob (j , observation, t,B, cell no) ) Fi [t, j] =arg_max(0,n-l, j , t,D, A) ,- i changes from 0 to n-1 endelse endfor endfor ,- termination for i=0, n-1 do begin

Z [i] =D [t_max-l, i] endfor ln_of_Prob_path=max(Z,k) _/ probability of the given path Q=intarr (t_max) _/ Q is the hidden state sequence

Q [t_max -l]=k for t=t_max-2 , 0, -1 do begin Q [t] =Fi [t+l , (Q lt + 1] ) ] endfor print , ' This is Q ' print , Q print , ' ' print , ' This is LogP of the hidden state sequence given the observation data : ' print , ln_of _Prob_path print , ' ' end

3 :begin print, 'button 3 pressed' for i=0, iteration_no do begin par_reestimation endfor print, 'this is new A' print, A print, ' ' print, 'This is new B' print, B print, ' ' print, 'this is new P' print, p print, ' ' end

6 : begin

_/ cell_no=2

,-w_s=3 ,-frame_no=15

,- t_max=frame_no - w_s + 1

,-yes_no_values=intarr (cell_no, frame_no)

,-yes_no_values [0 , *] =findgen (frame_no)

,-yes_no_values [1, *] =findgen (frame_no) _/ q=intarr (t_max)

,-q= [1,0, 0,1, 0,0, 1,0, 0,1, 0,0,1]

_/ correlation_state=l yes_no_values2=intarr (cell_no, frame_no) for i=0, frame_no -1 do begin for j=0, cell_no -1 do begin yes_no_values2 [j , i] =-1 endfor endfor for i=0, frame_no -1 do begin OK=0 for j=0, w_s -1 do begin /possibly change the j range to make the window smaller if(((i-j) ge 0) and ((i-j) It t_max) ) then begin if ( q(i-j) eq correlation_state ) then OK=l endif endfor for k=0, cell no -1 do begin if (OK eq 1) then Yes_no_values2 [k, i] =yes_no_values [k, i] endfor endfor

_/print, ' frmae number is', frame_no _/print, 't max is', t_max _/print, 'w_s is', w_s

xx=where (yes_no_values2 [0, *] eq -1, count) yes_no_values3=intarr (cell_no, frame_no - count) for j=0, cell_no -1 do begin k=0 for i=0, frame_no -1 do begin if (yes_no_values2 [j , i] ne -1) then begin yes_no_values3 [ j , ] = yes_no_values2 [ j , i] k=k+l endif endfor endfor

_/print, 'these are yes no values' _/print, yes_no_values

_/print, q

_/print, 'these are yes no values 2'

_/print, yes_no_values2

_/print, 'these are yes no values 3' _/print, yes_no_values3 print, 'State constrained yes and no values created' draw cross2 end endcase end pro par reestimation common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common markov, n,m, A, B,P, observation, f , g,d, Fi, Q, back, xsi, Gamma, n_i, a_i,b_i,p_i, c,w_s , t_max t_max=frame_no - w_s +1 _/ ! except=2

Let's recalculate forward variables: F f=dblarr (t_max, n) for i=0, n-1 do begin f (0, i) =P (i) *B (i, observation (0) ) / endfor for t=0, (t_max-2) do begin for j=0, (n-1) do begin _/ f (t, i) ,A(i, j ) ,-f (t+1, j ) =sum(0,n-l,t, j ,f ,A) *B (j , observation (t+1) ) ,- i goes from 0 to n-1

,- endfor

_/ endfor

,- 1 will apply scaling procedure to prevent underflows and increase accuracy

_/ c's will be the scaling coefficients c=dblarr (t_max) f_init=dblarr (t_max,n) f=dblarr (t_max,n) s=double (0) for i=0, n-1 do begin f_ini (0, i) =P (i) *B_prob (i, observation, 0,B, cell_no) s=s+f_init (0,i)

_/print, 'state is ', i, 'b_prob is ',

B_prob (i, observation, 0,B, cell_no) endfor c[0]=l/s for i=0, n-1 do begin f (0,i)=f_init (0,i) *c[0] endfor for t=0, (t_max-2) do begin s=double(0) for j=0, (n-1) do begin ,- f (t, i) , A(i, j ) f_init (t+1, j ) =sum(0,n-

1, t, j , f , A) *B_prob (j , observation, (t+1) ,B, cell_no) ,-i goes from 0 to n-

1 s=s+f_init (t+1, j ) endfor c[t+l]=l/s for j=0, (n-1) do begin f (t+1 , j ) =f_init (t+1 , j ) *c [t+1] endfor endfor

_/Let's define the backward variables: BACK

Back_init=dblarr (t_max,n)

Back=dblarr (t_max,n) for i=0, n-1 do begin

Back_init (t_max -l,i)=l

Back(t_max -1, i) =Back_init (t_max -1, i) *c [t_max-l] endfor for t=t_max-2, 0, -1 do begin for i=0, (n-1) do begin s=double (0) for j=0,n-l do begin s=s + A(i, j ) *B_prob (j , observation, (t+1) ,B, cell_no) *Back(t+l, j ^'. endfor Back_init (t, i) =s

Back(t,i) =Back_init (t,i) *c [t] endfor endfor

_/print, 'these are backward variables'

_/print, Back

_/print, 'space'

_/Now let's define the sequence of probabilities xsi(t,i,j) - it is the

,- robaility of being at state i at time t and at state j at time t+1 xsi=dblarr (t_max,n,n) for t=0,t_max-2 do begin s=double (0) for i=0,n-l do begin for j=0,n-l do begin s=s+F(t, i) *A(i, j) *B_prob (j , observation, (t+1) ,B,cell_no) *Back(t+l, j ) endfor endfor for i=0,n-l do begin for j=0,n-l do begin xsi [t,i, j] = (F(t,i) *A(i, j) *B_prob(j observation, (t+1) ,b,cell_no) *Back( t+l,j))/s endfor endfor endfor _/print, 'This is xsi'

_/print, xsi

_/print, 'space'

_/Gamma (t,i) is the probability of being at state i at t given the

_/observation and the model

Gamma=Dblarr (t_max, n) for t=0, t_max -1 do begin for i=0, n-1 do begin s=double(θ) for j=0, n-1 do begin s=s+xsi (t, i, j ) endfor

Gamma [t, i] =s endfor endfor

/print, 'This is gamma'

,-print , gamma

_/print, 'space'

/Now we will begin the reestimation procedures of initial parameters P,A,B

,-reestimation of P for i=0, n-1 do begin P (i) =gamma (1 , i) endfor

; Reestimation of A for i=0, n-1 do begin for j=0, n-1 do begin sl=double (0) s2=double (0) for t=0, t_max-2 do begin sl=sl+Xsi (t, i, j) s2=s2+Gamma (t, i) endfor

A(i, ) =sl/s2 endfor endfor _/Reestimation of B for j=0, n-1 do begin for k=0, m-1 do begin sl=double (0) s2=double (0) ,-for t=0, t_max-l do begin if (observation (t) eq k) then begin sl=sl+Gamma (t, j ) endif _. s2=s2+Gamma (t, j ) ,- endfor

B(j,k) =sl/s2 endfor endfor _/Reestimation of B for j=0, n-1 do begin for k=0, m-1 do begin sl=double (0) s2=double (0) for t=0, t_max-l do begin sl=sl + Gamma (t, j )* (observation [k, t] ) s2=s2+Gamma (t, j ) endfor

B(j,k)=sl/s2 endfor endfor

end

/ to sum n elements, used by hidden_markov_event . pro

function sum, x,y,t,j,f,A

_/s will be used to denote the sum s=f (t,θ)*A(θ, j) for i= (x+1) , y do begin s=s + f (t, i) *A(i, j) endfor return, s end APPENDIX C

_/DESCRIPTION: using this, the user finds the nonnegsvd by inputing the number of iterations and

_/the number of modes to keep. (Though it is called calculate_svd, it is a different linear _/factorization- nonnegative one.., but close to the general SVD) . _/MODIFICATION HISTORY:

_/ Created on April 23, 1999 by Vikram Sheel Kumar

_/ Last Modified on April 29, 1999 by Vikram Sheel Kumar pro calculate_svd common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common pp_com, tp_pixel_array base=WIDGET_BASE (/column, title=' Perform Non-negative SVD') draw=WIDGET_DRAW(base, xsize=400, ysize=400) modeskeep=CW_FIELD (base, title= 'Enter the number of modes that should be kept : ' , $ VALUE=cell_no, /integer, /return_events) niterations=CW_FIELD (base, title= 'Number of iterations:', $

VALUE=100, /integer, /return_events) button=WIDGET_BUTTON(base, value= ' Calculate SVD',- uvalue= ' button_pressed' ) widget_control, /realize, base state={modeskeep:modeskeep, niterations :niterations} widget_control, widget_info (base, /child), set_uvalue=state xmanager, ' calculate_svd' , base end _/ 1999, (cl Copyright by Vikram Sheel Kumar _/NAME: calculate_svd_event

_/DESCRIPTION: This is the event handler for calculate_svd.. It calls nonnegsvd_converge which calls nonnegfac_converge based on the _/number of iterations and the number of modes to keep. **Make sure to do all the required error checking...** _/MODIFICATION HISTORY: _/ Created on April 23, 1999 by Vikram Sheel Kumar

Modified on April 29, 1999 by Vikram Sheel Kumar _/ Last Modified on May 11, 1999 by Vikram Sheel Kumar to avoid possible conflict between this widget & visualization (as that needs ,-to have a slider built mkeep long, while mkeep can changes for different svd calculations) pro calculate_svd_event, event common pp_com, tp_pixel_array common svds_com, w, u, v, errors, mkeep, imu common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common flags, svd_calculated, cells_defined common vis_svdcom, visbase, visbase2 stateholder=WIDGET_INFO (event .handler, /Child) print, visbase, visbase2 vismanaged=WIDGET_INFO(visbase, /MANAGED) ,-if this is 1, the general visualization widget is open vismanaged2=WIDGET_INFO (visbase2, /MANAGED) _/if this is 1, the spatial visualization widget is open

WIDGET_CONTROL, stateholder, GET_UVALUE=state _/find what is in /Child

WIDGET_CONTROL, state .modeskeep, GET_VALUΞ=mkeep WIDGET_CONTROL, state . niterations , GET_VALUE=niter Widget_control, event. id, get_uvalue=uval

,-if ( (vismanaged eq 1) and (mkeep ne mkeepold) ) OR ( (vismanaged2 eq 1) and (mkeep ne mkeepold2))~ then begin

/number of modes being calculated is not number being showed in visualization.general / mess=WIDGΞT_MΞSSAGE( 'A visualization widget is open with a capacity for a different number of modes than you currently specify! Change the number of modes, or close the visualization widget and repeat . ' , /ERROR) if ( (vismanaged eq 1) OR (vismanaged2 eq 1) ) then begin mess=WIDGET_MESSAGE( 'To avoid conflicts with the number of modes specified, please shut your visualization widgets . ' , /INFORMATION) endif else begin if uval eq 'button_pressed' then begin if (mkeep gt cell_no) then begin mess=DIALOG_MESSAGE ( 'Number of modes cannot be more than the number of cells !', /ERROR) endif else begin nonnegsvd_converge, tp_pixel_array, w, u, v, mkeep, niter, errors imu=fltarr (x_size,y_size,mkeep) for i=0,mkeep-l do begin iml=fltarr (x_size,y_size) iml (location (*) . coord (0) , location (*) . coord (1) ) =u(*, i) _/reconstruct to dimensions of image for j=0, cell_no-l do begin imu (location (j ) . coord (0) - 3 : location (j ) . coord (0) +3 , location (j ) . coord (1) - 3 : location (j ) . coord (1) +3,i) =iml (location (j) . coord (0) , location (j ) .coor d(D) endfor endfor plot, w, xtitle='Mode Number', ytitle= ' Power

Contribution' svd_calculated=l endelse endif endelse end

_/ Name : cell_type_pt_svd.pro

,- Description: : This program calculates the svds for cells that have been previously selected by Multicell. It finds the points where the mouse _/ was clicked during Multicell and draws boxes of 3 pixels side around them .

; NOTE: This differs from the img_pt_svd programs as it allows you to specify the indices of the cells of a given type to perform an analysis only

,- on those specified cells (of the same type) . You can run this program separately for the different types of cells you see.

,- Variables: nx- size in x direction of images

_/ ny- size in y direction of images _/ location: saved Multicell file that has the coordinates of the chosen cells

,- ct : array of cell numbers of cells of the chosen type.

CELL 1 IS THE FIRST CELL ! !

_/ ima: array of pixel traces for the chosen cells ,- imu: spatial eigenmodes

_/ ts : temporal eigenmodes

_/ w: singular values

_/ ncomp: number of modes to look at. YOU CAN ONLY

HAVE AS^' MANY MODES AS NUMBER OF CELLS ANALYZED! ,- Modification History:

_/ Written by Vikram Sheel Kumar, Dec 20, '98

; Based on code written by Partha Mitra, Bell-Labs pro cell_type_pt_svd, nx, ny, ct, location, ima, imu, ts, w, ncomp ctype=ct-l _/as these must be indexed starting from

0. Input starts from cell 1 ncells=n_elements (ctype) nt=n_elements (ima (0, *) ) print, 'start subtraction of linear trend' xax=transpose (findgen (nt) ) nima=fltarr (ncells,nt) _/this will hold the calcium transients of only the cells specified by ctype for j=0,ncells-l do begin cell=ctype ( ) _/this is the cell of the co rrect type that has to be analyzed result=regress (xax, reform (ima (cell, *) ) , fltarr (nt) +l,yfit) ima (cell, *) =ima (cell, *) -yfit nima(j , *) =ima(cell, *) print, j endfor print, 'before svd' tqlisvd, nima, w, u, v, ncomp print, 'after svd' v=v (* , 0 : ncomp-1) u=u(*, 0 : ncomp-1) imu=fltarr(nx,ny, ncomp) for i=0,ncomp-l do begin iml=fltarr (nx,ny) iml (location (ctype) . coord (0) , location (ctype) . coord (1) ) =u(*, i)

/u has all the values for the spatial eigenmode. for j=0, ncells-1 do begin imu (location (ctype (j ) ) . coord (0) - 3 : location (ctype (j ) ) . coord (0) +3 , location (ctype (j) ) . coord (1) -

3 : location (ctype (j) ) . coord (1) +3,i) =iml (location (ctype (j) ) . coord (0) , lo cation (ctype (j ) ) . coord (1) ) endfor endfor ts=v end

NAME : choose_ce11s_svd

SYNOPSIS: choose_cells_svd

DESCRIPTION:

This is the program that creates the widget with which you can choose cells

_/ It asks you to open the frame you want to pick cells from and then uses the filepath for that image as the root directory to find the remaining images

This is called by svd_gui_event .

Created on April 16, 1999 (and deleted on April 16, 1999)

Recreated on April 23, 1999 by VSK

WIDGET TO CHOOSE

CELLS pro choose_cells_svd common mother_com, pixel_array, yes_no_values, coef, location, cell no, frame_no, file_name, direct,? total_frame_no, time_resolution, x_size, y_size, box_size common choose_cell_convert_delta, filepathl common choose_cells_com, locations_defined, variance_created locations_defined=θ variance_created=0 name=pickfile (/read, title='For Choosing Cells, Select Movie Slice', GET_PATH=filepathl) str_ing=name file_name=svd_stringname (str_ing) _/this is a procedure to extract the root from the full file name .. if (name eq ' ' ) then begin mess=WIDGET_MESSAGE ( 'No image was chosen', /INFORMATION) endif else begin

_/total_frame_no=find_frame_no (file_name, filepathl) _/this found total frames , but took TOO long first_frame=strcompress ( (file_name) + ( ' 000 ' ) ) image=tiff_read(name) _/this reads the first image to find the sizes size_o_i=size (image) x_size=size_o_i (1) _/this is the x size of the image y_size=size_o_i (2) ,-this is the y size of the image base3=WIDGET_BASE (/COLUMN, title='Cell Choosing Widget') draw5=WIDGET_DRAW (base3 , xsize=x_size, ysize=y_size) slide4=WIDGET_SLIDER (base3 , value=total__frame_no, maximum=total_frame_no, minimum=4, title='x: Include every other frame up to x to make variance array') slide2=WIDGET_SLIDER(base3, value=3 , maximum=60, minimum=l, title='y: Each box is 2y on a side') slide3=WIDGET_SLIDΞR (base3 , value=total_frame_no, maximum=total_frame_no, minimum=4, title= ' z : Calculate F/F over z frames ' ) buttonl=cw_bgroup (base3 , /row, ['Display Variance Array', 'Choose Cells', 'Find F/F'], /return_name, UVALUE=3)

WIDGΞT_CONTROL, /realize, base3 statel= {slide2 :slide2, slide:slide3 , slide4 : slide4, draw : drawδ }

WIDGET_CONTROL, WIDGET_INFO (base3 , /Child), SET_UVALUE=Statel xmanager, ' choose_cells_svd' , base3 endelse end

_/ SYNOPSIS: choose_cells_svd_event , event

_/ DESCRIPTION:

,- This is the event handler for the widget choose_cells_svd.pro

,- It allows the user to choose cells from a slice of a movie. It calls pick_box which is the program

,- all of the frames from the movie and first finds the pixel values and then converts them into F/F values. _/ MODIFICATION HISTORY

_/ Created on April 23, 1999 by VSK

_/ EVENT HANDLER OF THE

ABOVE pro choose_cells_svd_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common choose_cells_com, locations_defined, variance_created common choose_cell_convert_delta, filepathl common pp_com, tp_pixel_array common flags, svd_calculated, cells_defined common local, imvar

_/this makes it readable from 'Choose Cells' too stateholder=WIDGET_INFO (event .handler, /Child)

WIDGET_CONTROL, stateholder, GET_UVALUE=statel _/find what is in /Child cell_no=0B frame_no=0B if ( (size (event. alue) ) [1] eq 7) then begin _/ ,- i.e. if we've pressed a button which returns a string as its value... CASE event.value OF 'Display Variance Array' : BEGIN WIDGET_CONTROL, statel. draw, GET_VALUE=win

WIDGET_CONTROL, statel . slide4 , GΞT_VALUE=frame_no wset, win & cd, filepathl convert_alt_tiff , frame_no, file_name, imarr imav=2*total ( (imarr) ,3) /frame_no*l.00 imavsq=2*total ( (TEMPORARY (imarr^A2) ) , 3) /frame_no*l .00 imarr=0 imvar=imavsq-imav^A2 imavsq=0 & imav=0 poly=sfit (imvar, 2) ,-2D surface to normalize with imvar=imvar/poly & tvscl, imvar variance created=l

END

'Choose Cells' : BEGIN if (variance_created eq 0) then begin _/did not find variance image mess=WIDGET_MESSAGE('You must first display the variance array. ' , /ERROR) endif else begin tvscl, imvar ,-this clears if you want to rechoose cells WIDGΞT_CONTROL, statel. draw, GET_VALUE=win

WIDGET CONTROL, statel . slide2 , GET_VALUE=box_size location=pick_box(box_size, x_size, y_size) _/choose cells through this cell_no=location(0) .size _/ location (0) since we only want one size cells_defined=0 locations_defined=l mess=WIDGET_MΞSSAGE( 'Press Choose Cells to reselect cells. Otherwise press Find F/F', /INFORMATION) endelse

END

'Find F/F' : BEGIN if (locations_defined eq 0) then begin mess=WIDGET_MESSAGE( 'You have to choose cells before I can find the F/Fs!', /INFORMATION) endif else begin

WIDGET_CONTROL, statel . slide , GΞT_VALUE=frame__no cell_no=location(0) .size _/ location (0) since we only want one size pixel_array=convert_to_delta_f_box (location, frame_no, fi1e_name) tp_pixel_array=pixel_array WIDGET_CONTROL, event . top, /HOURGLASS cells_defined=l

_/this lets the programs know that cells have been defined! mess=WIDGET_MΞSSAGE( 'Calculated the F/Fs. NOW saving the locations of the cells...', /information) cell_location ,-this prints the cell locations to a file endelse END

ELSE: BEGIN WIDGET CONTROL, statel . slide2 , GET VALUE=box size

END ENDCASE endif end NAME : com

DESCRIPTION: this calculates the center of mass of a 2D array NOTE::im must be xsize*ysize*time pro com, im, xcm,ycm sz=size (im) xw=rebin (findgen (sz (1) ) , sz (1) , sz(2), /sample) yw=transpose (rebin (findgen (sz (2) ) , sz (2) , sz (1) , /sample) ) xcm=fltarr (sz (3) ) & ycm=xcm for i=0, sz(3)-l do xcm(i) =total (xw*im(*, *, i) ) /total (im(*, *, i) ) for i=0, sz(3)-l do ycm(i) =total (yw*im(*, *, i) ) /total (im(*, *, i) ) end

Name : convert_alt_tiff Synopsis: convert_tiff , rootname Description:

This opens many tiff files and converts them to one big array.

Modified by VSK on April 26, 1999 from convert_tiff_to_array so that it opens alternate frames pro convert_alt_tiff, frame_no, filename2, movie realname=strarr (4) realname (0) =filename2+string( ' 00 ' ) /if the frame_number is between 1 and 9 the units place of it's log will 0 and this will be called realname (1) =filename2+string ( '0 ' ) ,-if the frame_number is between 10 and 99, byte (it's log) =1 and this will be called realname (2) =filename2 ,-if the frame_number is between 100 and 999, byte (it's log) =2 and this will be called realname (3) =filename2 image=tiff_read (strcompress (filename2+string( ' 000 ' ) , /remove_all )) info=size (image) x=info(l) _/dimensions of image y=info(2) frames=floor (frame_no/2) print, 'too big' movie=bytarr (x, y, frames) movie (* , * , 0) =image for frame_number=2 , frame_no-l, 2 do begin image=tiff_read(strcompress (realname (byte (aloglO (frame_number) ) ) ÷string (frame_number) , /remove_all) ) movie (* , * , frame_number/2) =image image=0 endfor end

_/Name : convert_tiff_to_array

_/Synopsis: convert_tiff_to_array, rootname _/Description: This opens many tiff files and converts them to one big array.

Created by VSK on April 21, 1998 function convert_tiff_to_array, frame_no, filename2 realname=strarr (4 ) realname (0) =filename2+string (' 00 ' ) _/if the frame_number is between 1 and 9 / the units place of it's log will 0 and this will be called realname (1) =filename2+string( ' 0 ' ) ,-if the frame_number is between 10 and 99 ,

; byte (it's log) =1 and this will be called realname (2) =filename2 _/if the frame_number is between 100 and 999,

/ byte (it's log) =2 and this will be called realname (3 ) =filename2 image=tiff_read (strcompress (string (filename2) ÷string ( ' 000 ' ) , /remove_all) ) info=size (image) x=info(l) y=info (2) mov=bytarr (x,y, frame_no) mov(*, *, 0) =image ,- frame_number=l

,- image=tiff_read (strcompress (realname (byte (aloglO (frame_number) ) ) ÷string (frame_number) , /remove_all) ) ,- info=size (image) ,-x=info(l) _/dimensions of image

,-y=info (2)

,-mov=bytarr (frame_no, x, y) for frame_number=l, frame_no-l do begin print , frame_number image=tiff_read(strcompress (realname (byte (aloglO (frame_number) ) ) ÷stri ng(frame_number) , /remove_all) ) mov ( * , * , frame_number) =image endfor return, mov end

NAME: convert_to_de11a_f_box

SYNOPSIS : convert_to_delta_f_box, cell_location, frame_no, filename

DESCRIPTION:

This program converts the pixel values into deltaF/F values to normalize the data. ,- The deltaF/F value varies from 0 to 100%. By using this technique of measuring intensity, ,- a brighter cell will need a greater pixel change than a dimmer cell to achieve the same

_/ deltaF/F. The baseline is created by averaging the first 4 frames .

Formula used:

Delta F F - F(0)

X 100

F(0)

MODIFICATION HISTORY:

Created by Vikram Sheel Kumar in November, 1996

Modified on 12/3/96

Last Modified on 7/13/97 to complete documentation

function convert_to_delta_f_box, cell_location, frame_no, filename2 common choose__cell_convert_delta, filepathl common delta_f_main, baseline

!except=2 realname=strarr (4 ) realname (0) =filename2+string( ' 00 ' ) ,-if the frame_number is between _ and 9 the units place of it's log will 0 and this will be called realname (1) =filename2+string (' 0 ' ) _/if the frame_number is between 10 and 99 , byte (it's log) =1 and this will be called realname (2) =filename2 _/if the frame_number is between 100 and 999, byte (it's log) =2 and this will be called realname (3 ) =filename2 cell_no=cell_location(0) .size

/incorporate all this into a structure pixel_values=fltarr (cell_no, frame_no) delta_f_values=fltarr (cell no, frame_no) half_side=cell_location(0) .half_side

_/HALF LENGTH OF BOX SIDE total_side=half_side*2+l _/FOR EG. HALF_SIZE=3, TOTAL PIXELS CALCULATED OVER =7*7=49 area=total_side^A2

_/ OPEN EACH FRAME AND MEASURE THE PIXEL

VALUE OF EACH CELL THERE for frame_number=0, frame_no-l do begin tempo=tiff_read (strcompress (realname (byte (aloglO (frame_number) ) +string (frame_number) , /remove_all) ) for cell_numb=0, cell_no-l do begin x=cell_location (cell_numb) . coord (0) y=cell_location(cell_numb) . coord (1) box_array=tempo (x-half_side :x+half_side, y- half_side :y+half_side) pixel_values (cell numb, frame_number) = (total (box_array) ) /area _/this is the average pixel value in a box endfor endfor CONVERT EACH PIXEL ARRAY INTO AN ARRAY OF

DELTA F VALUES for cell_number=0 , cell_no-l do begin baseline=total (pixel_values (cell_number, 0:3))/4.00 _/average first 3 frames to get start delta_f_values (cell_number, *) = (pixel_values (cell_number, *) - baseline) / (baseline*l .0) endfor delta_f_values (*, *) =TΞMPORARY(delta_f_values (* , *))*100 return, delta_f_values end

,-Name : find_frame_no ^— _/Description: This procedure goes through a directory and finds the number of frames present. For example if the largest _/index it finds is Frame.343, it'll return 343 as the number of frames present function find_frame_no, file_name, filepath realname=strarr (4) realname (0) =file_name+strin (' 00 ' ) _/if the frame_number is between 1 and 9 the units place of it's log will 0 and this will be called realname (1) =file_name+string (' 0 ' ) _/if the frame_number is between

10 and 99, byte (it's log) =1 and this will be called realname (2) =file_name ,-if the frame_number is between

100 and 999, byte (it's log) =2 and this will be called realname (3) =file_name frame_number=0 name= (strcompress (realname (byte (aloglO (frame_number) ) ) ÷string (frame_n umber) , /remove_all) ) result=findfile (filepath+name, count=no_found) while (no_found gt 0) do begin frame_number=frame_number+l name= (strcompress (realname (byte (aloglO (frame_number) ) ) +string(f rame_number) , /remove_all) ) result=findfile (filepath+name, count=no_found) endwhile print , frame_number return, frame_number end pro img_svd, im, mask, imu, ts, w, ncomp nx=n_elements (im(*, 0, 0) ) ny=n_elements (im(0, *, 0) ) nt=n_elements (im ( 0 , 0 , *) ) iml=im(*,*,0)

,-mas =where (iml gt int_thr) nxl=n_elements (mask) ima=fltarr (nxl , nt) for i=0,nt-l do begin iml=im(*, *, i) ima (* , i) =iml (mask) endfor print, 'before svd' tqlisvd, ima, w, u, v, ncomp print, 'after svd' v=v(*, 0 :ncomp-l) u=u(*, 0 :ncomp-l) imu=fltarr (nx,ny, ncomp) for i=0, ncomp-1 do begin iml=fltarr (nx, ny) iml (mask) =u (* , i) imu(*, *, i) =iml endfor ts=v end

pro img_svd__rampsub, im, mask, imu, ts, w, ncomp nx=n_elements (im (* , 0 , 0) ) ny=n_elements (im(0, *, 0) ) nt=n elements (im(0, 0, *) ) iml=im(*, *, 0)

,-mask=where (iml gt int_thr) nxl=n_elements (mask) ima=fltarr (nxl,nt) for i=0,nt-l do begin iml=im(*, *, i) ima (*, i) =iml (mask) endfor print, 'start subtraction of linear trend' xax=transpose (findgen (nt) ) nmask=n_elements (mask) for j=θ,nmask-l do begin result=regress (xax, reform (ima (j , *) ) , fltarr (nt) +l,yfit) ima (j , *) =ima (j , *) -yfit print, j endfor print, 'before svd' tqlisvd, ima, w, u, v, ncomp print, 'after svd' v=v(*, 0 :ncomp-l) u=u(*, 0 :ncomp-l) imu=fltarr (nx, ny, ncomp) for i=0, ncomp-1 do begin iml=fltarr (nx, ny) iml (mask) =u(*,i) imu(*, *, i) =iml endfor ts=v end

_/Name : edsub . ro _/Description: This uses the running median filter to subtract the baseline _/ etc .

_/Variables: medsize is the size of the window used in calculating the filter

_/ cellsin is a cells X frames array that contains the unfiltered data

; cellsout has the cells after the base has been subtracted from them. Written by

Partha Mitra and Vikram Sheel Kumar, February 27 1999 Modified to take the absolute value instead of truncating. This is necessary for -ve spikes

pro medsub, cellsin, medsize, cellsout

cellsout=float (cellsin) ncells=n_elements (cellsin (*, 0) ) nframes=n_elements (cellsin (0, *) ) tmpl=fltarr (nframes+2*medsize) for i=0, ncells-1 do begin tmpl (medsize :medsize+nframes-l) =reform (cellsin (i, *) ) ,-tmpl= (tmpl-median (tmpl,medsize) ) >0 tmpl=abs (tmpl-median (tmpl,medsize) ) cellsout (i, *) =tmpl (medsize :medsize+nframes-l) endfor end pro nonnegfac, x^~, wl,vl

Iteration for computing nonnegative factorisation x = Matrix/ wl, vl = nonnegative actors, x = N x M wl = N x R vl = R x M wv=wl#vl xwv=x*0 nz=where (x gt 0) xwv(nz) =x(nz) /wv(nz) vl=vl* (transpose (wl) #xwv) xwv=xwv*0 wv=wl#vl xw (nz) =x (nz) /wv (nz) wl=wl* (xwv#transpose (vl) ) eps=le-20 norm=rebin (transpose (total (wl, 1) ) +eps,n_elements (wl (*, 0) ) ,n_elements ( wl(0,*))) wl=wl/norm end pro nonnegfac_converge, x wl,vl Iteration for computing nonnegative factorisation x = Matrix,- wl, vl = nonnegative factors, x = N x M wl = N x R vl = R x M wv=wl#vl xwv=x*0 nz=where (x gt 0) eps=le-20 xwv(nz) =x(nz) / (wv(nz) +eps) vl=vl* (transpose (wl) #xwv) xwv=xwv*0 wv=wl#vl xwv(nz) =x(nz) / (wv(nz) +eps) wl=wl* (xwv#transpose (vl) ) norm=rebin (transpose (total (wl, 1) ) +eps,n_elements (wl (*, 0) ) ,n_elements ( wl(0,*))) wl=wl/norm end pro nonnegsvd, mat, frac, u ^~, v^ nkeep, niter ,-frac is similar to the singular value matrix of before.. ,-nkeep is the no of modes to keep ,-niter= no of interations sz=size (mat) N=sz(l) M=sz (2) R=nkeep w=randomu (seed, N, R) v=randomu (seed, R, M) for i=0, niter do nonnegfac , ma , w, v recon=w#v error=fltarr (R) norm=total (mat^A2) for i=0,R-l do error(i) = total ( (mat-w(*, i) #v(i, *)) 2) /norm srt=sort (error) frac=l-error (srt) u=w(*, srt) v=v(srt, *) end

pro nonnegsvd_converge, mat, frac, u, v, nkeep, niter, errors _/svd, mat, w, u, v _/ index=reverse (sort (w) ) ,-w=w (index) _/u=u ( * , index) _;v=v(*, index) sz=size (mat) x=sz (1) y=sz (2) R=nkeep noise=randomu (seed,R,y) noise2=randomu(seed,x,R) w=fltarr (x, R) + .01*noise2 for i=0, R-l do w(i,i)=l v=mat (0 :R-1, *) +. l*noise

_;u=u(*,0:R-l)

;v=v(*,0:R-l)

_/for t=0, R-l do begin

,-u(*,t)=u(*,t)*sqrt(w(t) ) ,-v(*,t)=v(*,t)*sqrt(w(t) )

_/ endfor

,-w=abs (u)

,-v=abs (transpose (v) ) errors=fltarr (niter+1) for j=0, niter do begin nonnegfac_converge , mat , w, v recon=w#v error=fltarr (R) norm=total (mat^A2) for i=θ,R-l do error(i) = total ( (mat-w(*, i) #v(i, *)) ^A2) /norm srt=sort (error) frac=l-error (srt) u=w(*, srt) v=v(srt,*) errors (j ) =total ( (mat-recon) 2) /norm endfor end

_/Name : pick_box Version4.0 _/Synopsis: picks_box, cell_numbers

/Description: Through this function, the user picks the cells he/she wants . These coordinates can then ,-be input into the get_cells_oplot_win function which gives the pixel values and binary transformation _/of the cells.

7

/UPDATED VERSION AS OF 6/25/96 by Vikram Sheel Kumar /this program is a version of pick_cells that highlights a cell as soon as a point in the cell is _/clicked.

7 _/picks_box can only be done after pick_cells has been done so the boxes can be _/chosen from the same image as the cells were chosen

7 _/OLD statements that might be imp. later: ,- naming the files : -

NOTE, CAN'T READ FILES WITH LAST NAME 'TIFF' names=realname ( 0 ) +strcompress ( (string (format= ' (i) ' , f indgen (f ram e number) +1) ) , /remove_all )

7 names=filename (0) +strcompress ( (string (format= ' (i) ', findgen (fram e_number) +1) ) , /remove_all)

,- names=filename (0) +strcompress ( (string (format= ' (i, " .tiff ") ' , findgen (frame_number) +1) ) , /remove_all) _/ data being saved: - _/ file_name2=string ( ("Saved_Locations_Box@") , strcompress (strmid (systime () , 3, 10) , $

/remove_all) , '_' , strcompress (strmid (systime () , 14,2) ) ) close, 3 openw, 3, FILEPATH(file_name2, subdirectory=dir) _/OPEN LUN TO WRITE INTO

,- printf, 3, 'half side of box=', box_size

/ printf, 3, 'For cell ', i+l, 'you clicked at { ' , cell (i) .coord,

'}'

_/ close, 3 _/CLOSE THE LUN

MAIN PROGRAM

BEGINS function pick_box, box_size, xsize, ysize

celll=replicate ( {struc, coord: intarr (2) } , 500) ,-Stru to hold the location (x and y) of MAX NO OF CELLS CHOOSABLE=500

_/ for the first cell boxlength=box_size*2 cell_number=l cursor, x, y, /up, /device celll (cell_number) . coord (0)=x celll (cell_number) . coord (1) =y

_; give an error if the mouse was not clicked on the image if ( (x ge xsize-boxlength OR x le boxlength) OR (y ge ysize- boxlength OR y le boxlength) ) then begin mess=WIDGET_MΞSSAGE( 'You clicked outside the allowed area! Type xmanager to return to Multicell.', /Error) retail endif

,- to make sure the first cell doesn't fail the while test we make cell(O) diff from cell(l) celll(O) .coord(0)=-l for every other cell while (x ne celll (cell_number-l) . coord (0) OR y ne celll (cell_number-l) .coord (1) ) do begin

,-if these cond fail, you have double clicked on a cell /and the program stops counting cells.

_; draw the box and number for j=-box_size, box_size do begin plots, [x+j, x+j], [y-box_size, y+box_size] ,

/device _/plots lets me draw points endfor xyouts, x-boxlength*3 , y+boxlength, cell_number, $ _/xyouts lets you write on an image alignment=.5, /device cell_number=cell_number+l

. read the cursor position from the screen cursor, x, y, /up, /device

. give an error if the mouse was not clicked on the image if ( (x ge xsize-boxlength OR x le boxlength) OR (y ge ysize- boxlength OR y le boxlength) ) then begin mess=WIDGΞT_MESSAGE( 'You clicked outside the allowed area! Type xmanager to return to Multicell.', /Error) retail endif

. p_U the x and y coordinates into the structure celll (cell_number) . coord (0) =x celll (cell_number) .coord (1) =y endwhile make another structure with the correct no of cells i.e cell number-1 cell=replicate ( {struct, coord : intarr (2) , size:0, half_side: 0.00} , cell_number-l) for i=0, cell_number-2 do begin cell(i) . coord=celll (i+l) .coord endfor cell . size=cell_number-l _/total number of cells is cell_number-l since one extra cell is cell .half_side=box_size _/recorded when you double click to stop counting, return, [cell] end

_/Name: preprocess_wid

,- Decription: This is the widget that allows you to preprocess your data (smooth, medsub, abs) /Modification History

/ Created on 4/16/99 by Vikram Sheel Kumar to add to svd_gui pro preprocess_wid common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size base3=WIDGET_BASE (/column, title= ' Preprocess ' ) draw5=WIDGET_DRAW(base3,xsize=550, ysize=300) slide2=WIDGET_SLIDΞR (base3 , minimum=l , maximum=cell_no, title='Cell Number') button=cw_bgroup (base3 , /row, [ 'Plot ', 'Median Subtract', 'Smooth', 'Take Absolute Value¹, 'Return Original Values ' ] , /return_name, UVALUΞ=0) medwin=CW_FIELD (base3 , title= 'Running Window Width for Median Subtraction' ,VALUE=20, /INTEGER, /RETURN_EVENTS) smoothwin=CW_FIELD (base3 , title=' Width of Smoothing Window ',VALUE=3, /INTEGER, / ETURN_ΞVENTS) state={slide:slide2 , medwin:medwin, smoothwin:smoothwin, draw : draws } widget_control, base3 , /realize widget_control, base3 , set_uvalue=state xmanager, 'preprocess_wid' , base3 end

Name: preprocess__wid_event Description: Widget handler Modification History

Created on 4/16/99 by Vikram Sheel Kumar pro preprocess_wid_event, event common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no , time_resolution, x_size, y_size, box_size common pp_com, tp_pixel_array

Widget_control, event. top, get_uvalue=state3 WIDGET_CONTROL, state3. slide, GET_VALUE=cell_numb WIDGET_CONTROL, state3. smoothwin, GET_VALUE=swin WIDGET_CONTROL, state3.medwin, GET_VALUE=mwin WIDGET_CONTROL, state3.draw, GET_VALUE=win_num if (mwin le 1) then begin mess=WIDGET_MESSAGE( 'Median filter window width must be greater than 1 ! ' , /ERROR) endif else begin if ( (swin It 2) OR (swin ge frame_no) ) then begin print , frame_no mess=WIDGET_MESSAGE ( 'Smoothing window width must be >2 and <Total number of frames ! ' , /ERROR) endif else begin if ( (size (event.value) ) [1] eq 7) then begin _{/ /} i.e. if we've pressed a button which returns a string as its value... Case event.value Of

'Smooth' : BEGIN wset, win_num for i=0, cell_no-l do begin tp_pixel_array ( i , * ) =smooth (tp_pixel_array ( i , * ) , swin) endfor plot , f indgen ( frame_no) *time_re solution, tp_pixel_array (cell_numb- l , * ) END

' Median Subtract ' : BEGIN wset , win_num medsub , tp_pixel_array, win, tp_pixel_array plot , f indgen (frame_no) *time_resolution, tp_pixel_array (cell_numb-l , * ) END

'Take Absolute Value' : BEGIN wset, win_num tp_pixel_array=abs (tp_pixel_array) plot, findgen (frame_no) *time_resolution, tp_pixel_array (cell_numb-l, *) END

'Return Original Values': BEGIN wset, win num tp_pixel_array=pixel_array plot, findgen(frame_no) *time_resolution, tp_pixel_array (cell_numb-l, *) END 'Plot' : BEGIN wset, win_num plot, findgen (frame_no) *time_resolution, tp_pixel_array (cell_numb-l, *)

END

ELSE: ENDCASE endif else begin wset, win num plot, findgen (frame_no) *time_resolution, tp_pixel_array (cell_numb-l, *) endelse endelse endelse end

_/ Name : pt_img_svd_rampsub

,- Description: : This program calculates the svds for cells that have been previously selected by Multicell. It finds the points where the mouse

_/ was clicked during Multicell and draws boxes of 3 pixels side around them.

Variables: nx- size in x direction of images ny- size in y direction of images _/ location: saved Multicell file that has the coordinates of the chosen cells

_/ ima: array of pixel traces for the chosen cells imu: spatial eigenmodes ts : temporal eigenmodes w: singular values ncomp: number of modes to look at Modification History:

Written by Vikram Sheel Kumar, Nov 24, '98 Based on code written by Partha Mitra, Bell-Labs pro pt_img_svd_rampsub, nx, ny, location, ima, imu, ts, w, ncomp ncells=n_elements (ima (*, 0) ) nt=n_elements (ima ( 0 , *) ) print, 'start subtraction of linear trend' xax=transpose (findgen (nt) ) nmask=n_elements (mask) for j=0,ncells-l do begin result=regress (xax, reform (ima (j , *) ) , fltarr (nt) +l,yfit) ima (j , *) =ima (j , *) -yfit print, j endfor print, 'before svd' tqlisvd, ima, w, u, v, ncomp print, 'after svd' v=v(*, 0 :ncomp-l) u=u(*, 0 :ncomp-l) imu=fltarr (nx,ny, ncomp) for i=0, ncomp-1 do begin iml=fltarr (nx,ny) iml (location (*) . coord (0) , location (*) . coord (1) ) =u(*, i) for j=0, ncells-1 do begin imu (location (j ) . coord (0) - 3 : location (j ) . coord (0) +3 , location (j ) . coord (1) -

3 : location (j ) . coord (1) +3 , i) =iml (location (j ) . coord (0) , location (j ) . coor d(l)) endfor endfor ts=v end

SYNOPSIS: r_noise, seed, length DESCRIPTION: This creates random noise of length 'length' to be included in a transient of our test data. MODIFICATION HISTORY:

Created on January 4, 1999 by Vikram Sheel Kumar function r noise, length rl=randomu(seed, length) plot, rl r2=float (randomu(seed, length) *length) plot, r2 rr=rl(r2) _/this scrambles the random array return, rr end Name : rampsub

Synopsis: rampsub, image

Description: This program subtracts linear ramps from movies.

Variables: image=the 3-D array that holds the frames (x,y) and time as the z coordinate _/ Created by Partha Mitra function rampsub, im ,- subtracts linear ramps sz=size (im) nx=sz (1) _/the x length of frames ny=sz(2) _/the y height of frames nt=sz(3) ,-the time or frame number imout=bytarr (nx,ny,nt) /make a floating type array of the same size (this will take memory!) for i=0,nx-l do begin ,-go through every point on the image print, i for j=θ,ny-l do begin xax=transpose (findgen (nt) ) ,-transose this for the regression equation result=regress (xax, reform(im(i, j , *) ) , fltarr (nt) +l,yfit) imout (i, j , *) =im(i, j , *) -yfit endfor endfor return, imout end

,- Name: randomize

,- Synopsis:

,- Description: This procedure takes the pixel traces of the chosen cells and shifts each

,- of them by a random amount. We do this so we can compare the svds for a random set and real data set .

/ Variables :

_/ ima: array of pixel traces for the chosen cells _/ it is ncells X nframes in dimensions

,- Modification History

,- Created by Vikram Sheel Kumar, Nov 24, 98

,- Modified by VSK on Dec 14, 98 to improve the randomizing function by scrambling the random nos / Last Modified by VSK on March 23, 1999 to change the name to

'randomize' function randomize, ima nframes=n_elements (ima(0,*) ) ncells=n_elements (ima (*, 0) ) rima=ima rarra=randomu (seed, ncells) _/make an array with ncell random entries rin=floor (randomu (seed, ncells) *ncells) _/make another array to randomly subscr ipt the prev. one rar=rarra(rin)

_/ scramble the array for i=0, ncells-1 do begin rima (i, *) =shift (ima (i, *) , rar (i) *nframes) endfor return, rima end

pro read_cret, fname, nframe, dat ,-dat is the output array

_/this is supposed to read the photodiode data with a carriage return dat=fltarr (nframe) close, 1 openr, 1, fname tmp= ' ' readf , 1 , tmp dat (0) =tmp for i=l, nframe-1 do begin readf, 1, tmp ,-to get rid of the carriage return readf , 1 , tmp if (tmp eq -3374.23) then print, i da (i) =float (tmp) endfor close, 1 end function read3col, f ame, ncell, ntim dat=fltarr (ncell , ntim) tmp= ' ' openr, 1 , fname for i=0,ncell-l do begin for j=0,ntim-l do begin readf , 1 , tmp wherek=strpos (tmp, 'K' ) dat (i, j ) =float (strtrim(strmid (tmp, wherek+1, 8) ) ) endfor endfor close, 1 return, dat end function readdata, fname , ncell, ntim dat=fltarr (ncell , ntim) ,- tmp= ' ' openr, 1, fname for i=0,ncell-l do begin _/print, ' cell ' , i for j=0,ntim-l do begin readf , 1 , tmp dat (i, j) =tmp endfor endfor close, 1 return, dat end

Name: convert_tiff_to_array

Synopsis : convert_tiff_to_array, rootname

Description:

This opens many tiff files and converts them to one big array.

-, Created by VSK on April 21, 1998 function rebin_tiff_to_array, frame_no, filename2 realname=strarr (4) realname (0) =filename2+string( ' 00 ' ) ,-if the frame_number is between 1 and 9 the units place of it's log will 0 and this will be called realname (1) =filename2+string (' 0 ' ) ,-if the frame_number is between 10 and 99, byte (it's log) =1 and this will be called realname (2) =filename2 ,-if the frame number is between 100 and 999, byte (it's log) =2 and this will be called realname (3) =filename2 image=tiff_read (strcompress (filename2+string( ' 000 ' ) , /remove_all )) info=size (image) x=info(l) _/dimensions of image y=info(2) movie=fltarr (x/2, y/2, frame_no) image=tiff_read (strcompress (filename2+string ( ' 000 ' ) , /remove_all) ) movie (*,*, 0) =rebin(image, x/2, y/2) for frame_number=l, frame_no-l do begin image=tiff_read(strcompress (realname (byte (aloglO (frame_number) ) ) ÷string (frame_number) , /remove_all) ) movie (*, *, frame_number) =rebin (image, x/2, y/2) endfor return, movie end

7+

_/Name : remove_drift .pro _/Synopsis: remove_drift, imarr _/Introduction: This is an automation of the drift removing techniques . ,- It finds the drift and then removes it via interpolation. _/ Description of Variables : imarr is a 3D array with

_/ imarr (frame_number, x, y) . It holds the frames of the movie we are analyzing. Modification History

Created by Vikram Sheel Kumar on Jul 23, 98

pro remove_drift, imarr size_info=size (imarr) if (size_info (0) ) It 3 then

MΞSSAGE_WINDOW ( ' The image array must be 3D¹, /ERROR) endif else begin frame_no=size_info (1) the number of frames in movie xsize=size_info (2) the x size ysize=size_info (3) the y size imagetl=bytarr (xsize, ysize) ,-this holds the image (t) imaget2=bytarr (xsize, ysize) ,-this holds image (t+1) for t=0, frame_no-2 do begin _/inspect the pairs of frames imagetl (* , *) =imarr (t, * , *) ,- image (t) imaget2 (*, *) =imarr (t+1, *, *) _/image (t+1) imagetl=rebi (imagetl, 2*xsize, 2*ysize) _/pad up 4 imaget2=rebin(imatet2 , 2*xsize, 2*ysize) /pad up 4 ,-perform ffts on images fftl=fft (imagetl) fft2=fft(imaget2) fft2=conj (fft2) correl_f=fftl*fft2 /this is the correlation fnt peak_loc=where (correl_f eq max(correl_f) ) print, loc endfor end

pro run_lme0, y, n, y_lme dn=l nt=n_elements (y) nwin=fix( (nt-n) /dn+1) tim=findgen (n) w=replicate (l,n) yfit=fltarr (nwin, n) for j=0,nwin-l do begin tmp=regress (transpose (tim) , reform (y(dn*j :dn*j+n-l) ) , w, tmpfit) yfit (j , *) =reform(tmpfit) endfor y_line=fltarr (nt) norm=fltarr (nt) for j=0,nwin-l do begin for t=0,n-l do begin indx=j *dn+t y_line (indx) =y_line (indx) +yfit (j ,t) norm(indx) =norm(indx) +1.0 endfor endfor mask=where (norm gt 0) y_line (mask) =y_line (mask) /norm(mask) end pro str_mgname common test , str_ing, state3 x=strpos (str_ing, ' : ' , 0) while (x ne -1) do begin x=x+l str_ing=strmid(str_ing, (x) ) x=strpos (str_ing, ' : ' , 0) endwhile x=strpos (str_ing, ' . ' , 0) if (x ne -1) then begin str__ing=strmid(str_ing, 0,x+l) endif

,-help, str_ing end

Name : superimposed_plots

Synopsis :

Description: This program plots the deltaF/F values of a single cell per time.

Modification History:

Created by Vikram Sheel Kumar on October 23, 1996 Last Modified by Vikram Sheel Kumar on May 4, to incorporate into / svd_gui and change name to ' superimposed_plots ' pro superimposed_jplots, delta_f_values, cell_number, frame_number, time_resolution

!p.multi=0 window, 2

,-plot, findgen (frame_number) *time_resolution, delta_f_values (0, *) , yrange=[-50, 100+ (cell_number*20) ] , $ title='All cells', ytitle= ' F/F¹, xtitle= ' time in seconds ' plot, findgen (frame_number) *time_resolution, delta_f_values (0, *) , yrange= [-20,20] , $ title='All cells', ytitle= ' F/F', xtitle='time in seconds ' xyouts, -6, 0, 1 if cell_number gt 1 then begin for cell_numb=l, cell_number-l do begin

/ oplot, findgen (frame_number) *time_resolution, (delta_f_values (cell_numb, *) + (cell_numb*20) ) oplot, findgen (frame_number) *time_resolution, (delta_f_values (cell_numb, *) ) /xyouts, -6, cell_numb*20, cell_numb+l a= (delta_f_values (cell_numb, *) eq max ( (delta_f_values (cell_numb, *) ) a=where (a eq 1) xyouts, a,max( (delta_f_values (cell_numb, *) ) ) , cell_numb+l endfor endif end

NAME : SVD_gui

SYNOPSIS: SVD_gui

DESCRIPTION: This widget will be used to load data, pre-process, perform an SVD, and then visualize the modes created. pro svd_gui common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common flags, svd_calculated, cells_defined common vis svdcom, visbase, visbase2 visbase=long(-l) visbase2=long(-l) svd_calculated=0 set flags to 0 cells_defined=0 base=WIDGET_BASE (/column, title='SVD GUI') stuff= { cw_pdmenu_s , flags :0, name:¹'} details=[{ cw_pdmenu_s, 0, 'Previous Experiment' cw_pdmenu_s, 0, 'Choose Cells' }, $ cw_pdmenu_s , 0, ¹ Superimpose Traces ' } , $ cw_pdmenu_s , 0, 'Preprocess ' } , $ cw_pdmenu_s , 0, ' SVD ' } , $ cw_pdmenu_s , 1, 'Visualize' } , $ cw_pdmenu_s , o, ' General ' } , $ cw_pdmenu_ s , 2, 'With Locations'}] pull_down=cw_pdmenu (base, details, /return_full_name, UVALUE=12) stuff2= {cw_pdmenu_s2, flags :0, name:'¹} details=[{ cw_pdmenu_s, 0, 'Save Analysis'}, $

{ cw_pdmenu_s , 0 , ' Color ' } , $

{ cw_pdmenu_s, 0, 'Retail'}, $

{ cw_pdmenu_s , 0 , ' Exit ' } ] pull_down2=cw_pdmenu(base, details, /return_full_name, UVALUE=13]

frameid=CW_FIELD (base, title= ' otal Number of Frames ' ,VALUE=total_frame_no, /INTEGER, $

/RETURN_EVENTS) timeid=CW_FIELD (base, title='Time Resolution in seconds' ,VALUE=time_resolution, /FLOATING, /RETURN_EVENTS) state= { frame :frameid, time :timeid} widget_control , /realize, base widget_control, widget_info (base, /Child), SET_UVALUE=state xmanager, ' svd_gui ' , base end

Name : svd_gui_event ^'

Synopsis: svd_gui_event , event

Description: This is the event handler for the svd_gui window. Modification History:

Created by Vikram Sheel Kumar on April 21, 99 pro svd_gui_event , event common trtother_com, pixel_array, yes_no_values , coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common pp_com, tp_pixel_array common flags, svd_calculated, cells_defined stateholder=WIDGET_INFO (event .handler, /Child) WIDGET_CONTROL, stateholder, GET_UVALUE=state ,-find what is in /Child

WIDGET_CONTROL, state. frame, GET_VALUΞ=total_frame_no WIDGET CONTROL, state. time, GET_VALUE=time_resolution Case event.value OF

' Previous Experiment ' : BEGIN tempstate=state data_file=dialog_pickfile (/read, title= 'Select the saved variables file', GET_PATH=filepathl) if (data_file eq ' ' ) then begin mess=WIDGET_MESSAGΞ ( 'No saved variables file specified', /INFORMATION) endif else begin restore, data_file mess=WIDGET_MESSAGE( 'Done loading variables',

/INFORMATION)

WIDGET_CONTROL, tempstate . frame, SET_VALUΞ=total_frame_no

WIDGET_CONTROL, tempstate . time , SET_VALUE=time_resolution endelse END

' Choose Cells ' : BEGIN if ( (total_frame_no eq 0) or (time_resolution eq 0) ) then begin mess=WIDGET_MESSAGΞ ( 'You must enter the total number of frames and time resolution before proceeding!', /ERROR) endif else begin choose_cells_svd endelse END

¹ Superimpose Traces ' : BEGIN if ( (total_frame_no eq 0) or (time_resolution eq 0)) then begin mess=WIDGET_MESSAGE( 'You must enter the total number of frames and time resolution before proceeding!', /ERROR) endif else begin if (cells_defined eq 1) then begin superimposed_plots, tp_pixel_array, cell_no, frame_no, time_resolution endif else begin mess=WIDGET_MΞSSAGE( 'Choose cells before plotting.', /ERROR) endelse endelse END

'Preprocess': BEGIN if ( (total_frame_no eq 0) or (time_resolution eq 0)) then begin mess=WIDGET_MESSAGE('You must enter the total number of frames and time resolution before proceeding!', /ERROR) endif else begin if (cells_defined eq 1) then begin preprocess_wid endif else begin mess=WIDGET_MESSAGE(' Choose cells before plotting. ' , /ERROR) endelse endelse

END

' SVD ' : BEGIN if ( (total_frame_no eq 0) or (time_resolution eq 0) ) then begin mess=WIDGET_MESSAGE ( 'You must enter the total number of frames and time resolution before proceeding!¹, /ERROR) endif else begin if (cells_defined eq 1) then begin calculate_svd endif else begin mess=WIDGET_MESSAGEp Choose cells before calculating the SVD.', /ERROR) endelse endelse

END

'Visualize. General ' : BEGIN if ( (total_frame_no eq 0) or (time_resolution eq 0) ) then begin mess=WIDGET_MESSAGE ( 'You must enter the total number of frames and time resolution before proceeding!', /ERROR) endif else begin if (svd_calculated eq 1) then begin visualize_wid endif else begin mess=WIDGET_MESSAGE ( 'You must calculate the SVD before visualizing the modes .', /ERROR) endelse endelse

END

'Visualize.With Locations': BEGIN if ( (total_frame_no eq 0) or (time_resolu ion eq 0)) then begin mess=WIDGET_MESSAGE( 'You must enter the total number of frames and time resolution before proceeding!', /ERROR) endif else begin if (svd_calculated eq 1) then begin visualize_wid_locs endif else begin mess=WIDGET_MΞSSAGE( 'You must calculate the SVD before visualizing the modes .', /ERROR) endelse endelse END

' Color' : xloadct

'Retail ' : retail

'Save Analysis': BEGIN data_file=pickfile (/write, file= 'Variables_svd#_' , title=' Create the saved variables file') if (data_file eq ' ' ) then begin mess=WΣDGET_MESSAGE( 'No saved variables file specified', /INFORMATION) endif else begin

SAVE, /VARIABLES, FILENAME=data_file, all, /verbose mess=WIDGET_MESSAGE( 'Data from this experiment saved', /INFORMATION) endelse END

'Exit' :BEGIN

WIDGΞT_CONTROL, /DESTROY, event . top END ENDCASE end

_/ Name: svd_stringname _/ Synopsis: svd_stringname (str_ing)

_/ Description: This finds the root of a file from its long name. It looks for the part of the file before .XYX where XYZ goes from 000 to 9999 etc.. Inputs: str_ing is the full name of the file Modification History:

Written by Paul Vartan, Summer 98 _/ Last Modified by Vikram Sheel Kumar, April 29, 99 to make it a function rather than a procedure

function svd_stringname, str_ing x=strpos (str_ing, ' : ' , 0) while (x ne -1) do begin x=x+l str_ing=strmid(str_ing, (x) ) x=strpos (str_ing, ' : ' , 0) endwhile x=strpos (str_ing, ' . ' , 0) if (x ne -1) then begin str_ing=strmid(str_ing, 0,x+l) endif return, str_ing end Copyright (c) 1996-1997 by Vikram Sheel Kumar

NAME: svd_to_delta_f_box

SYNOPSIS: svd_to_delta_f_box, cell_location, frame_no, filename

DESCRIPTION:

This program converts the pixel values into deltaF/F values to normalize the data.

,- The deltaF/F value varies from 0 to 100%. By using this technique of measuring intensity,

,- a brighter cell will need a greater pixel change than a dimmer cell to achieve the same

,- deltaF/F. The baseline is created by averaging the first 4 frames .

Formula used: Delta F F(0)

X 100

F(0)

MODIFICATION HISTORY:

Created by Vikram Sheel Kumar in November, 1996 Modified on 12/3/96

Modified on 7/13/97 to complete documentation Last Modified on 5/4/99 to move to svd_gui and start calculating from .001 instead of .000 function svd_to_delta_f_box, cell_location, frame no, filename2 common choose_cell_convert_delta, filepathl common delta_f_main, baseline

!except=2 realname=strarr (4) realname (0) =filename2+string( ' 00 ' ) ,-if the frame_number is between 1 and 9 the units place of it's log will 0 and this will be called realname (1) =filename2+string( ' 0 ' ) _/if the frame_number is between

10 and 99, byte (it's log) =1 and this will be called realname (2) =filename2 _/if the frame_number is between

100 and 999, byte (it's log) =2 and this will be called realname (3) =filename2 cell no=cell location (0) . size

; incorporate all this into a structure pixel_values=fltarr (cell_no, frame_no) delta_f_values=fltarr (cell no, frame_no) half_side=cell_location(0) .half_side

_/HALF LENGTH OF BOX SIDE total_side=half_side*2+l _/FOR EG.

HALF_SIZE=3, TOTAL PIXELS CALCULATED OVER =7*7=49 area=total side 2 OPEN EACH FRAME AND MEASURE THE PIXEL

VALUE OF EACH CELL THERE for frame_number=l, frame_no do begin tempo=tiff_read(strcompress (realname (byte (aloglO (frame_number) ) ) ÷string (frame_number) , /remove_all) ) for cell_numb=0, cell_no-l do begin x=cell_location (cell_numb) . coord (0) y=cell_location(cell_numb) . coord (1) box_array=tempo (x-half_side :x+half_side, y- half_side :y+half_side) pixel_values (cell_numb, frame_number) = (total (box_array) ) /area ,-this is the average pixel value in a box endfor endfor CONVERT EACH PIXEL ARRAY INTO AN ARRAY OF

DELTA F VALUES for cell_number=0, cell_no-l do begin baseline=total (pixel_values (cell_number, 0:3))/4.00 _/average first 3 frames to get start delta_f_values (cell_number, *) = (pixel_values (cell_number, *) baseline) / (baseline*l .0) endfor delta_f_values (*, *) =TΞMPORARY(delta_f_values (*, *) ) *100 return, delta_f_values end pro tqlisvd, mt, w, u, v, num

Description: calculates the SVD of a skinny matrix (M=pxq, p>>q) usual IDL svd is slower in this case (I think due to bad IDL programming - not sure why) . Variables: mt- matrix

m=mt nx=n_elements (m (* , 0) ) ny=n_elements (m(0, *) ) if ny gt nx then m=transpose (m) nu=n_elements (m(*, 0) ) nel=n_elements (m ( 0 , * ) ) mr=fltarr (nel, nel) for i=0,nel-2 do begin for j=i+l,nel-l do mr (i,j ) =total (float (m(*, i) ) *float (m(*,j ) endfor mr=mr+transpose (mr) for i=θ,nel-l do mr (i, i) =total (float (m(*, i) ) *float (m(*, i) ) ) ;m=0 tred2 ,mr,d, e tqli, d,e,mr ds=reverse (sort (d) ) d=abs (d) d=sqrt(d(ds) ) v=fltarr (nel , num) for i=0,num-l do v(*, i) =mr (*, ds (i) ) u=fltarr (nu,num) for i=θ,num-l do u (*, i) =m#v(*, i) /d(i) w=d if ny gt nx then begin tmp=u u=v v=tmp endif end

,- Description: This is the widget that shows the spatial and temporal modes of the data. _/ Modification History

_/ Created on April 29, 1999 by Vikram Sheel Kumar pro visualize_wid common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common svds_com, w, u, v, errors, mkeep, imu common wid, num_plots common vis_svdcom, visbase, visbase2 num_plots=θ ,-no plots have been done, so the first must be a plot, not Oplot visbase=WIDGET_BASE (/column, title= 'Visualization Widget') drawl=WIDGET_DRA (visbase, xsize=500 ,ysize=100) draw2=WIDGET_DRAW (visbase, xsize=500,ysize=150) slide=WIDGET_SLIDER (visbase, value=l, maximum=mkeep , minimum=l, $ title='Mode Number', uvalue= ' slider_pressed' ) buttonl=WIDGET_BUTTON (visbase , value= ' Plot ' , uvalue= ' plot_pressed ' ) draw3=WIDGET_DRAW (visbase, xsize=500,ysize=150) slide2=WIDGET_SLIDER (visbase, value=l, maximum=cell_no, minimum=l, $ title='Cell Number', uvalue=' slider2_pressed' ) button2=WIDGET_BUTTON (visbase, alue= 'Add Cell to Plot ' , uvalue= ' Add ' ) button3=WIDGΞT_BUTTON (visbase , alue= ' Clear Cells ' , uvalue= ' Clear ' )

WIDGET_CONTROL, /realize, visbase loadct, 27 statel={drawl: rawl, draw2:draw2, draw3:draw3, slide: slide, slide2 :slide2, button2 :button2 , $ button3 :button3 } WIDGET_CONTROL, visbase, set_uvalue=statel xmanager, 'visualize_wid' , visbase end

,- Created on April 29, 1999 by Vikram Sheel Kumar pro visualize_wid_event, event common svds_com, w, u, v, errors, mkeep, imu common wid, num_plots common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common pp_com, tp_pixel_array

WIDGET_CONTROL, event. top, get_uvalue=state WIDGET_CONTROL, state . drawl , get_value=winl WIDGΞT_CONTROL, state. draw2, get_value=win2 WIDGET_CONTROL, state . draw3 , get_value=win3 WIDGET_CONTROL, state. slide, get_value=mode_no WIDGET_CONTROL, state . slide2 , get_value=plotcell WIDGET_CONTROL, event. id, get_uvalue=uval !y.minor=-l _/this is so that you can see the values on the y axis- suppress the tick-marks if ( (string (uval) eq 'plot_pressed' ) OR (string (uval) eq ' sliderjpressed' ) ) then begin wset, win2 & plot, v(mode_no-l, *) , xtitle= 'Frames ' , ytitle= ' Amplitude ' , title= ' Temporal Mode ' wset, winl & plot, findgen (cell_no) +1, u (*,mode_no-l) , psym=2, xtitle='Cell Number', title=' Spatial Mode¹ num_plots=0 endif if (string (uval) eq 'Add') then begin _/Add Cell to Plot was pressed if (num_plots eq 0) then begin wset, win3 & plot, tp_pixel_array (plotcell-1, *) , title= ' Superimposed Plots ' num_plots=l endif else begin wset, win3 & oplot, tp_pixel_array (plotcell-1, *) , color=num_plots*100 num_plots=num_plots+l endelse endif if (string(uval) eq 'Clear') then begin _/Clear Cells was pressed num_plots=0 wset, win3 & erase endif !y. minor=0 end

Description: This shows the spatial modes of the components with specific locations ,- Modification History:

_/ Created on April 29, 1999 by Vikram Sheel Kumar _/ Last Modified on May 11, 1999 by VSK to include mkeepold2 to avoid conflicts with the svd wid pro visualize_wid_locs common mother_com, pixel_array, yes_no_values, coef, location, cell_no, frame_no, file_name, direct, $ total_frame_no, time_resolution, x_size, y_size, box_size common svds_com, w, u, v, errors, mkeep, imu common vis_svdcom, visbase, visbase2 visbase2=WIDGET_BASE (/column, title= ' Scroll through the spatial modes ' ) draw=WIDGET_DRAW(visbase2, xsize=x_size, ysize=y_size) ' slide=WIDGET_SLIDER(visbase2, value=l, maximum=mkee , minimum=l, title='Mode Number', $ uvalue= ' slider_pressed' ) buttonl=WIDGET_BUTTON (visbase2 , value= ' Plot ' , uvalue= ' plot_jpressed ' )

WIDGET_CONTROL, /realize, visbase2 statel= {draw: draw, slide: slide} WIDGET__CONTROL, visbase2, set_uvalue=statel xmanager, 'visualize_wid_locs ' , visbase2 end

Name: visualize_wid_locs_event

Description: This will show the spatial locations of the cells from the modes ; Modification History: _/ Created on April 29, 1999 by Vikram Sheel Kumar pro visualize_wid_locs_event, event common svds_com, w, u, v, errors, mkeep, imu

WIDGET_CONTROL, event. top, get_uvalue=state WIDGET_CONTROL, state. draw, get_value=win WIDGET_CONTROL, state. slide, get_value=mode_no WIDGET_CONTROL, event. id, get_uvalue=uval if ( (string (uval) eq 'plot_pressed' ) OR (string (uval) eq ' sliderjpressed¹ ) ) then begin loadct, 5 wset, win & tvscl, imu (*, *,mode_no-l) endif end

Name : zerodeny_2d

Synopsis : zeroderiv_2d, m

Description of Variables: m is the matrix over which we want to find the singular points

Introduction: This function uses regress to find the coefficients of the polynomial fit to matrix. Then it takes the first derivative to find the local min and max. For the unlikely event of the det being zero, you will get an arithmetic error.

Modification History

Written by P. Mitra on 7/28/98 function zeroderiv_2d, m

,- Compute position of maximum in matrix m _/ use biquadratic interpolation sz=size (m) nx=sz(l) ,-nx is the number of columns in the matrix ny=sz(2) ,-ny is the number of rows in the matrix x=reform(rebin(findgen (nx) ,nx,ny, /sample) ,nx*ny)

,-we do interpolation via rebinning using the nearest neighbor

,-then we reform it into one long column y=reform (rebin (transpose (findgen (ny) ) ,nx,ny, /sample) ,nx*ny) ml=reform(m,nx*ny) _/reform this into a long column xx=fltarr (5,nx*ny) /for the regression, we want 5 coeff and have nx*ny pts . xx(0,*)=x & xx(l,*)=y & xx(2,*)=x 2 & xx(3,*)=y 2 & xx(4,*)=x*y a=regress (xx, ml, fltarr (nx*ny) +1.0) /the last term is the weight=l for all

,- finding the stationary pts .

,-we have to find the first derivatives w.r.t x and y _/and simultaneously solve them ,-a (0) +2a(2)x+a(4)y=0 ,-a(l)+a(4)x+2a(3)y=0 det=4*a(2) *a(3) -a (4) 2 xmax=(a(4)*a(l) -2*a(3) *a(0) ) /det ymax=(a(4)*a(0) -2*a(2) *a(l) ) /det return, [xmax, ymax] end

Claims

1. A method for analyzing a sequence of data arrays comprising the steps of:

(a) selecting at least one type of region of interest of said data arrays; (b) selecting at least one region of interest for each of said at least one selected type of region of interest selected in step (a) from said data arrays;

(c) transforming said sequence of data arrays into a first simplified data array having a first dimension equal to the number of said at least one region of interest selected in step (b), a second dimension equal to the number of data arrays in said sequence, and a third dimension equal to the number of said selected types of regions of interest;

(d) detecting events of interest in said at least one selected region of interest;

(e) storing said detected events of interest in a second simplified data array of binary data, having the same dimensions as said first simplified data array; and

(f) analyzing data in one of the data arrays selected from the group consisting of said first simplified data array and said second simplified data array, to determine relationships between said detected events of interest.

2. The method of claim 1, wherein said sequence of data arrays comprises a series of two dimensional arrays corresponding to a time series of two dimensional observations of an object.

3. The method of claim 1, wherein said object is neuronal cell slice containing a plurality of neurons and wherein said step of selecting a region of interest comprises selecting one or more neurons in said series of two dimensional observations.

4. The method of claim 2, wherein one type of region of interest is selected in step (a) and wherein said step of selecting at least one region of interest comprises selecting at least one portion of said sequence of data arrays that corresponds to a particular spatial region of said two dimensional observations of said object.

5. The method of claim 4, wherein said transforming step comprises extracting data from each array in said sequence of data arrays that corresponds to said particular spatial region of interest and storing said data in said first simplified array, each data entry having an associated first dimensional index and an associated second dimensional index, wherein said first dimensional index is between one and the number of spatial regions corresponding to selected selected regions of interest of said sequence of data arrays, and said second dimensional index is between one and the number of two dimensional observations of said obj ect in said sequence of data arrays.

6. The method of claim 5, wherein said step of detecting events of interest comprises calculating a statistical mean and statistical standard deviation from a data population consisting of all entries in said first simplified array having identical first dimensional indexes, for each of said first dimensional indexes; determining for each entry in said first simplified array whether said entry exceeds, by a predetermined number of said standard deviation associated with said entry, the mean associated with said entry and denominating such a data entry an event.

7. The method of claim 6, wherein entries in said second simplified data array have the same associated first and second dimensional indexes as corresponding entries in said first simplified data array and wherein said storing said detected events of interests comprises storing a one in said second simplified array when the data entry with the corresponding first and second dimensional indexes in said first simplified array is denominated an event and storing a zero in said second simplified array otherwise.

8. The method of claim 5 , wherein said detecting events of interest comprises deteπnining whether a first data entry in said first simplified array exceeds, by a threshold amount, a second data entry in said first simplified array wherein said second data entry has an identical first dimensional index as said first data entry and a second dimensional index corresponding to an earlier point in time than said first data entry and denominating said second data entry an event.

9. The method of claim 8, wherein entries in said second simplified data array have the same associated first and second dimensional indexes as corresponding entries in said first simplified data array and wherein said storing said detected events of interests comprises storing a one in said second simplified array when the data entry with the corresponding first and second dimensional indexes in said first simplified array is denominated an event and storing a zero in said second simplified array otherwise.

10. The method of claim 1 , wherein said step of analyzing comprises plotting at least a portion of said data in said first simplified data for visual analysis.

11. The method of claim 1 , wherein said step of analyzing comprises detecting said events of interest that are coactive and determining whether the number of coactive events is statistically significant.

12. The method of claim 11, wherein said step of detecting events of interest that are coactive comprises detecting instances where said events of interest are detected in two or more of said regions of interest simultaneously.

13. The method of claim 11 , wherein said step of detecting events of interest that are coactive comprises detecting instances were events of interest are detected in two of said regions of interest simultaneously at a plurality of locations along said second dimension of said second simplified data array.

14. The method of claim 1, wherein said step of analyzing comprises calculating a strength of correlation between at least two regions of interest based on the number of coactive events of interest occurring in said at least two regions of interest and displaying a correlation map illustrating the strength of correlation between said regions of interest by lines connecting said regions of interest wherein the thickness of each of the lines is proportional to said calculated strength of correlation between respective regions of interest connected by the line.

15. The method of claim 1 , wherein said step of analyzing comprises displaying a cross-correlogram between events of interest occurring in at least one of said regions of interest.

16. The method of claim 1 , wherein said step of analyzing comprises detecting at least one hidden Markov state sequence from said second simplified data array.

17. The method of claim 16, wherein said step of analyzing further comprises displaying a cross-correlogram between events of interest occurring in one of said regions of interest while said region of interest is in one of said detected hidden Markov states.

18. The method of claim 1, further comprising: before said step of detecting, performing a singular valued decomposition on said data in said first simplified data array to calculate a predetermined number of eigenmodes; modifying said data in said first simplified data array by removing the data that corresponds to the first of said predetermined number of eigenmodes; and storing said modified data into said first simplified data array.

19. A method for analyzing a sequence of data arrays comprising the steps of:

(a) selecting at least one type of region of interest of said data arrays;

(b) selecting at least one region of interest for each of said at least one type of region of interest selected in step (b) from said data arrays;

(c) transforming said sequence of data arrays into a first simplified data array having a first dimension equal to the number of said selected regions of interest, a second dimension equal to the number of data arrays in said sequence, and a third dimension equal to the number of said selected types of regions of interest; and

(d) performing a singular valued decomposition on said first simplified data array to determine relationships between said regions of interest.

20. A system for analyzing a sequence of data arrays comprising in a data processor: a region type selector for selecting at least one type of region of interest of said data arrays; a region selector for selecting at least one region of interest for each of said selected type of region of interest from said data arrays; a first data transformer for transforming said sequence of data arrays into a first simplified data array having a first dimension equal to the number of said selected regions of interest, a second dimension equal to the number of data arrays in said sequence, and a third dimension equal to the number of said selected types of regions of interest; an event detector for detecting events of interest in said regions of interest; a second data transformer for storing said detected events of interest into a second simplified data array of binary data, having the same dimensions as said first simplified data array; and a data analyzer for analyzing data in one of the data arrays selected from the group consisting of said first simplified data array and said second simplified data array, to determine relationships between said detected events of interest.

21. The system of claim 20, wherein said sequence of data arrays comprises a series of two dimensional arrays corresponding to a time series of two dimensional observations of an object.

22. The system of claim 21, wherein said object is neuronal cell slice containing a plurality of neurons and wherein said step of selecting a region of interest comprises selecting one or more neurons in said series of two dimensional observations.

23. The system of claim 21, wherein said region type selector selects one type of region of interest and wherein said region selector selects at least one portion of said sequence of data arrays that corresponds to a particular spatial region of said two dimensional observations of said object.

24. The system of claim 23, wherein said region selector receives coordinates corresponding to said particular spatial region from a user.

25. The system of claim 23 , wherein said first data transformer extracts data from each array in said sequence of data arrays that corresponds to said particular spatial region of interest and stores said data in said first simplified array, wherein each data entry in said first simplified array has an associated first dimensional index and an associated second dimensional index, wherein said first dimensional index is between one and the number of spatial regions corresponding to selected regions of interest of said seeuence of data arrays and said second dimensional index is between one and the number of two dimensional observations of said object in said sequence of data arrays.

26. The system of claim 25, wherein said event detector furthe comprises: a statistical calculator for calculating a statistical mean and statistical standard deviation from a data population consisting of all entries in said first simplified data array having identical first dimensional indexes, for each of said first dimensional indexes; and a comparator for determining for each entry in said first simplified array whether said entry exceeds, by a predetermined number of said standard deviations associated with said entry the mean associated with said entry and denominating such a data entry an event.

27. The system of claim 26, wherein said second data transformer stores entries in said second simplified data array having the same associated first and second dimensional index as corresponding entries in said first simplified data array and wherein said second data transformer stores a one in said second simplified array when the data entry with corresponding first and second dimensional indexes in said first simplified array is denominated an event and stores a zero in said second simplified data array otherwise.

28. The system of claim 25, wherein said event detector determines whether a first data entry in said first simplified array exceeds, by a threshold amount, a second data entry in said first simplified data array, wherein said second data entry has an identical first dimensional index as said first data entry and a second dimensional index corresponding to an earlier point in time than said first data entry, and denominates said second data entry an event.

29. The system of claim 28, wherein said second data transformer stores entries in said second simplified data array having the same associated first and second dimensional index as corresponding entries in said first simplified data array and wherein said second data transformer stores a one in said second simplified array when the data entry with corresponding first and second dimensional indexes in said first simplified array is denominated an event and stores a zero in said second simplified data array otherwise.

30. The system of claim 20, wherein said data analyzer plots at least a portion of said data in said first simplified array for visual analysis.

31. The system of claim 20, wherein said data analyzer detects said events of interest that are coactive and determines whether the number of coactive events is statistically significant.

32. The system of claim 31 , wherein said data analyzer detects said events of interest that are coactive comprises operability to detect instances where said events of interest are detected in two or more of said regions of interest simultaneously.

33. The system of claim 31 , wherein said data analyzer detects instances where events of interest are detected in two of said regions of interest simultaneously at a plurality of locations along said second dimension of said second simplified data array.

34. The system of claim 20, wherein said data analyzer calculates a strength of correlation between at least two regions of interest based on the number of coactive events of interest occurring in said at least two regions of interest and displays a correlation map illustrating the strength of correlation between said regions of interest by lines connecting said regions of interest wherein the thickness of each of the lines is proportional to said calculated strength of correlation between respective regions of interest connected by the line.

35. The system of claim 20, wherein said data analyzer displays a cross- correlogram between events of interest occurring in at least one of said regions of interest.

36. The system of claim 20, wherein said data analyzer detects at least one hidden Markov state sequence from said second simplified data array.

37. The system of claim 36, wherein said data analyzer displays a cross- correlogram between events of interest occurring in one of said regions of interest while said region of interest is in one of said detected hidden Markov states.

38. The system of claim 20, wherein said system further comprises: a decomposer to perform a singular valued decomposition on said data in said first simplified data array and calculating a predetermined number of eigenmodes; a data modifier for modifying said data in said first simplified data array by removing the data that corresponds to the first of said predeteπnined number of eigenmodes; and data storage for storing said modified data into said first simplified data array.

39. A system for analyzing a sequence of data arrays comprising in a data processor: a region type selector for selecting at least one type of region of interest of said data arrays; a region selector for selecting at least one region of interest for each of said selected type of region of interest from said data arrays; a first data transformer for tranforming said sequence of data arrays into a first simplified data array having a first dimension equal to the number of said selected regions of interest, a second dimension equal to the number of data arrays in said sequence, and a third dimension equal to the number of said selected types of regions of interest; a decomposer for performing a singular value decomposition on said first simplified data array to determine relationship between said regions of interest.