Let's think about a new culture of science.

Academic production before internet

Today, science is produced by the following basic infrastructure. This simplistic picture depicts how humans in a lab interact with the nodes making up the infrastructure. Humans basically use computational resources to analyze data by writing code. The input to these nodes overall is extremely sparse. That is we generally do not have other peoples' data available, we generally cook our own code specifically for our own data, we generally use our own computational resources.

Unit infrastructure for scientific production

This unit is the fundament for academia. Research is carried out by the same infrastructure that is simply replicated across geographies and time. Of course, there could be labs collaborating with each other, of course we could be using an external grid engine to run our tasks, of course we might download a toolbox to run analyses. These are all connections that are not shown in this picture. 

But the emphasis here is that we spent most of our time to reinvent the wheel by

-writing the same piece of code that many people had done it in the past, 

-collecting yet a new dataset instead of generating a new hypothesis compatible with available datasets, 

-buying large computers that could be used by other people

-hiring system administrators that are doing exactly the same work as in another lab

the list is long...

Each lab in this culture becomes a specialized idiosyncratic creature with its own way of doing things. Politically this implies committing in long-term fixed-costs to maintain an academic infrastructure that is short-sighted and benefits mainly to the labs short-term agenda.  Academia mainly benefits from the contributions of labs in the form of publications, which is considered as the unique currency in academic reward system. What is the impact of this system on the society? In the light of current replication crises in science, it is hard to be optimistic.

Infrastructure for academic work. Culture of mine.

This type of infrastructure organization has mainly historical reasons. This model is archaic, and has been a good model for the pre-internet era where people and systems were connected sparsely with each other, where it made sense to travel to a conference and meet other people.

A new way of doing science at the age of cloud-based systems

We have to rethink about how to place boxes shown in the previous pictures, how to set novel incentive mechanisms, and how to organize the work flow across scientists and nodes. Let's talk about this simple picture.

A novel infrastructure for academic work. Culture of sharing.

Outsourcing the storage and compute resources to a cloud service (e.g. AWS, GDC or some supranational public cloud service yet to be put in place) are for the benefit of the society in terms of reducing overall costs. 

However the main point here is not about outsourcing storage and compute resources. The real reason for this move is for making datasets accessible to other scientists. And in the long-term this actually means making data to be publicly available to all citizens.

The only thing that is specific to a given lab is the data that is collected there. That's what labs should do: collect data. Most importantly data must to be stored according to strict standardization. That is to each data set, a map has to be associated, that will help people on how to navigate this data set. Furthermore, every dataset should be stored with a minimal code that ensure basic access to data. Also, most often datasets spans multiple modalities. For example, my fMRI datasets are typically bundled together with pupil recordings and heart-beat recordings. We therefore need not only a standardization for storage of specific kind of datasets, but we also need a way to create dataset-bundles that represents an experiment in a flexible manner. A principled way to bundle standardized datasets. Let's call this step 1.

The other thing that labs do is to write code to process their data. To my opinion this is where the biggest challenge is located, namely on finding a system where people can collaborate and create something together. Assuming that the step 1 is solved, the code that is written will also be publicly available. Therefore, code that is written will be directly connected to a dataset type. 

For example, if I am trying to detect peaks in a more or less periodical physiological recording, I will not start looking for literature, find someone's algorithm, implement my version of it. I will simply search for code that is compatible with this type of data, browse among alternative codes, read comments to figure out strengths and weaknesses, consider ratings and incorporate that code to my pipeline.

Basically putting up an analysis will be about creating a pipeline using previously coded nodes or coding new nodes when the analysis has not been previously carried out. When something doesn't work as expected, code needs to be improved via collaboration. Writing good quality code will be one great novel incentive for scientists.

Another challenge is to find a way how to fund this novel system. This is certainly beyond the capacity of a single start-up. This is also beyond the scope of a single lab or institute. I also don't think today's national states are visionary enough to take such moves. To my opinion this could only be established by some tech giants who have the know-how required to solve all these problems.

The culture of "mine" in science at the age of cloud-based analysis systems*

Scientists needs computing power and storage space for their data sets. For scientific institutions, this translates onto long-term fixed-costs that are relatively high. Resources required for buying hardware, keeping network infrastructure, paying system administrators to take care of these masses of electronics is a considerable overhead. As a result, public scientific institutes spend lots of money and human resources to create and maintain infrastructures for storing and analyzing scientific data sets.

Yet, one scientific institute is pretty much the exact replica of another one when it comes to hardware demands. That is, resources that are needed in one place should in theory be very similar to another place. Therefore, instead of investing money for system administrators, storage and computational resources, scientific institutes may actually lease these services from cloud-based infrastructures with more flexible pricing opportunities and lack of overhead. Replacing your system administrator with two PhD students is an appealing idea after all.

There is actually nothing illuminating in this view because this has been actually happening already since more than 10 years in the corporate world. Many hosting companies offers as also VNC based system to connect to their servers and use software on powerful machines. Beyond simple hosting companies, Google Cloud Computing and Amazon AWS making the transformation real by integrating all sort of compute, storage, parallelization tools and selling it as a service.

Where are we in neuroscience? Some important milestones are becoming finally a reality in natural sciences, I think that the point of no return is also being slowly reached for neuroscience. I believe this because standardization procedures on how to store and share datasets is becoming more and more mainstream, and this shift has the potential to change day-to-day scientific enterprise radically. For example, Open Neuro is one such platform, where you can upload your brain imaging dataset using the BIDS format, and let analyses run on these servers. I think this is just start of a big scale transformation on how we do science.

Here is how I think how:

(1) Scientific publication

The way we publish our reports didn't change probably since the times of Fisher or even Newton. The world today is a very different place, but many of the novel tools that have been invented in the internet-based communication era have not been incorporated into the way how we conduct science today. OK, instead of sending a manuscript to the editor's office via post, we are today using emails, fine.

For example, the scientific reviewing system did not incorporate crowd-sourcing mechanisms to evaluate the quality of scientific papers. The decision of whether a manuscript or research proposal is worth being published stays largely within the hands of few not-randomly selected referees and an editor. The process is opaque, prone to biases and has no means to stop formation of small-world cartels that mutually benefit from positive biases.

The re-distribution of reputation is not based on metrics that reflect the long-term value of person for science in general. In the best case, reputation is equivalent to your h-index, which is heavily biased by the random success of your publication track, not how good a scientist your are. Metrics that ensures long-term advancements of science are typically not included. For example, we lack a metric that judges a professor based on the number of students that became also professors in the last 5-10 years. The infrastructure to achieve a better and more democratic system is in place since more than a decade. I believe this change will come faster with cloud-based systems decreasing the cost for storage and computational resources.

In the very near future, I believe any serious publication will also need to contain the related datasets, the analysis pipeline and make it publicly available to all scientific community (but also other citizens). This is already happening, and many journals let you agree with their terms of sharing data promptly when requested. However, the definition of "prompt" is also very subjective. For example, you may want to read this twit-storm to see a recent example. Even if it was obligatory to upload the dataset, the re-evaluation of the data is not within the responsibilities of the referees. This means that modifying an existing system incrementally to make it more and more suitable for the current demands of scientific democratization is not enough, we need a radically different way of publishing science.

When the data is stored and analysis ran in a cloud-based system, there will be no more excuses for reviewers for not being involved in the data analysis, as the time it will take for them to have a closer look on the data and the analysis pipelines will be insignificant. Therefore, I believe that any serious publications will take the concerted efforts of, on the one hand authors who designed the experiment, collected data and wrote the initial draft of the paper, and on the other, reviewers who will be required to contribute in the data analysis using infra-structures provided by the cloud-based storage and computational infrastructures. There will possibly be not much difference between collaborators of today and reviewers of tomorrow.

(2) Cloud-based analysis

Most of published reports use similar methods, which are re-invented again and again by generations of PhD and postdoc crew, which is a complete waste of time and resources. I believe actually there could possibly not be a more inefficient system than today's science. A large company would not be able to function like this.

Once we start talking about cloud-based storage and analysis pipelines, it will also be possible to run these analyses automatically on a server. You will need to tick the checkbox for this or that analysis and receive the results as an email in the form of a presentation or a web page (example) to click/browse around. This is of course an over simplification, but what I would like to say is that scientists will spend more time on (1) standardizing their datasets to be able to run analyses on the cloud-based system and (2) making analysis pipelines that are compatible with standardized datasets. Therefore, many scientists will use this time to record more data.

(3) End of culture of "mine"

One of the most intriguing anthropological traits of the daily scientific enterprise, is what I call the culture of "mine". This is not something that is somewhere out there, it is right inside our offices. By this I mean the way how students, PhDs, postdocs and professors (the whole crew basically) are closed to the idea of sharing and opening their projects to external influences. Most often if not always, a project is assigned to a single person in the lab, and this person is expected to run this project until the end. Because the person believes that it is her/his project, he/she can control the monopoly together with his/her boss on how this project has to run and adjust the level of external factors (politics). This results in a very conservative set of interactions between people, as any request of help, or any communication can be seen as a contribution to the project. The culture of mine, will of course be there and start the appropriate set of behaviors to not let this happen. Unfortunately, there are countless examples of authorship disputes which appear exactly from this type of culture.

Once the opportunity to upload your dataset and run your analysis in a cloud-based system is within the reach, there will be no reason to not open your data and let other people analyze it in ways different than what you have actually thought would be most appropriate. In a crow-sourced science, you will own your data, but will actually allow other people to look into it. Pretty much the same way, when people are allowed to look at you when you are walking in the street. The constructive discussions that follows during this process belongs to all parties and can be moderated by the person who created the dataset. I believe there will be a shift in the way how people conceptualize the way how they own projects and data, replacing culture of mine with crowd-sourced intelligence.

I found this article from Jeremy Freeman, entitled "Open source tools for large-scale neuroscience" which made me super happy as it expresses many of the thoughts I scratched on this post in a systematic and professional manner.

*This article has a bias from the perspective of a neuroscientist.

Standardization of fMRI analysis with SPM using the FancyCarp Matlab toolbox

Fancycarp is a repository for streamlining fMRI preprocessing and analysis based on SPM using Matlab.

With the FancyCarp toolbox, you create a folder hierarchy for your fMRI project. Download the fMRI data from the DICOM server (only for internals of IFSN). Preprocess fMRI data with SPM. Conduct first- and second-level analyses.

It can be integrated together with SCR, Pupil anaylses (experimental) for the analysis of autonomic responses.

How do I proceed?

1/ Get the repo

First clone this repository and switch to mrt/main branch.

onat@neocortex:/tmp$ git clone https://github.com/selimonat/fancycarp.git
Cloning into 'fancycarp'...
remote: Counting objects: 2154, done.
remote: Compressing objects: 100% (4/4), done.
remote: Total 2154 (delta 0), reused 0 (delta 0), pack-reused 2150
Receiving objects: 100% (2154/2154), 986.53 KiB | 801.00 KiB/s, done.
Resolving deltas: 100% (1414/1414), done.
Checking connectivity... done.

onat@neocortex:/tmp/fancycarp$ git checkout mrt/main
Branch mrt/main set up to track remote branch mrt/main from origin.
Switched to a new branch 'mrt/main'

If you would like to version-track your repository (which you should) create a new branch following /mrt/XXXX, where XXXX is the name of your project.

2/ Add paths

Fire up Matlab, add Fancycarp to Matlab's path. Don't forget to add SPM to your repository.

addpath('/tmp/fancycarp')

3/ Project Object

Fancycarp is based on OOP. This is convenient as it helps one to define a Project with a set of properties and methods, and encapsulate them in one single file. Functions that are specific to a project are declared in the ProjectObject, whereas other functions, which are for example specific to participants are declared accordingly at their respective objects. This removes clutter and helps defining clear functional separation of code. ProjectObject contains all sort of functions used for low-level house keeping. For example, finding the path to a subject at a certain run you can run the following code.

>> Project().pathfinder(2,2)
ans =
/tmp/mynextsciencepaper/data//sub002/run002/

The role of the ProjectObject is to provide functionality that spans multiple participants at once. For example, second-level analysis is defined in the ProjectObject, as this is a project-wide operation that uses all participants. Similarly, path finding, DICOM transformations, sanity checks, Mysql requests are all defined here at the ProjectObject.

Now, when you start a new project the first thing is to adapt the properties defined in the Project object to your project. These are defined in the file called Project.m.

properties (Hidden, Constant)%adapt these properties for your project
        %All these properties MUST BE CORRECT and adapted to one owns project
        ====> Path to your project folder.
        path_project          = '/tmp/mynextsciencepaper/data/';
        ====> Path to your SPM installation.
        path_spm              = '/home/onat/Documents/Code/Matlab/spm12-6685/';
        ====> Here enter your participants number.
        trio_sessions         = {'PRISMA_19873' 'PRISMA_19875'};
        ====> For every participant enter the number of acquisition runs of TRIO/PRISMA scanner to be downloaded.
        dicom_serie_selector  = {[8 19 20 21 6 7 17 18 ] [8 19 20 21 6 7 17 18 ] };
        ====> Runs [8, 19, 20...] will be distributed to folders [1, 2, 3, ...]
        dicom2run             = repmat({[1:8]},1,length(Project.dicom_serie_selector));
        ====> For field map corrections, enter their runs as stored in the dicom server.
        runs_fieldmap         = [{5 6} {7 8}];%The order is important: first one is the magnitude and the second one is the phase.
        ====> Related to VDM correction. 
        apply_vdm             = [{1}   {2 3 4}];        
        data_folders          = {'midlevel' 'mrt' 'design'};%if you need another folder, do it here.
        ====> Your TR, High-pass filtering values, wheterh to exclude surface extraction (requires CAT12)
        TR                    = 0.99;              
        HParam                = 128;%parameter for high-pass filtering
        surface_wanted        = 0;%do you want CAT12 toolbox to generate surfaces during segmentation (0/1)                
        smoothing_factor      = 4;%how many mm images should be smoothened when calling the SmoothVolume method        
    end

4/ Data Folders

Standardization of data folders in a project is the first step to facilitate code sharing and reproduction. To create a folder hierarchy, run

>> Project().CreateFolderHierarchy()

This will use Properties of the ProjectObject and loop over all participants and runs to create a folder structure to download anatomical and functional data. Check it out with the shell command tree.

The community seems to be settled upon the BIDS (brain imaging data structure) standard. BIDS will also be adopted in future versions of Fancycarp toolbox.

>> !tree                      
.
└── data
    ├── sub001
    │   ├── run000
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run001
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run002
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run003
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run004
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run005
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run006
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run007
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   └── run008
    │       ├── design
    │       ├── midlevel
    │       └── mrt
    ├── sub002
    │   ├── run000
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run001
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run002
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run003
    │   │   ├── design
    │   │   ├── midlevel
    │   │   └── mrt
    │   ├── run004
    │   │   ├── design

5/ Downloading Anatomical Data (Internal)

Now we have the folder hierarcy, we can proceed with downloading the anatomical and functional data from the internal server. To download data of a given participant, we need to create a Subject instance. Thanks to OOP inheritance, SubjectObject receives all the methods defined in the ProjectObject in addition to the methods that itself defines. Running the Subject().get_hr method will download the subject's anatomical scans, convert that to NifTi and store them at thesubXXX/run000/mrt folder.

>> s = Subject(1)
Subject Constructor for id:1 is called:
s = 
  Subject with properties:
              id: 1
            path: '/tmp/mynextsciencepaper/data//sub001/'
    trio_session: 'PRISMA_19873'
       total_run: 8
>> s.get_hr                                 
get_hr:
Will now dump the latest HR (16:05:01)
Dicom Server returns:
=====
Database: prisma
Patient: XXXXXX
  * Examination: PRISMA_XXXXX (XXXX, F)                          [ 1|  23|  5192]
    + Study:  1 (2017-11-08 11:15:00)                            [  |  23|  5192] 
      - Series:   1 {localizer (expectb-a, blank)              } [  |    |    11] 
      - Series:   2 {AAHead_Scout_64ch-head-coil               } [  |    |   128] 
      - Series:   3 {AAHead_Scout_64ch-head-coil_MPR_sag       } [  |    |     5] 
      - Series:   4 {AAHead_Scout_64ch-head-coil_MPR_cor       } [  |    |     3] 
      - Series:   5 {AAHead_Scout_64ch-head-coil_MPR_tra       } [  |    |     3] 
      - Series:   6 {gre_field_map, 2mm, filter M              } [  |    |    90] 
      - Series:   7 {gre_field_map, 2mm, filter M              } [  |    |    45] 
      - Series:   8 {ep2d_bold, mb3, loc                       } [  |    |  1193] 
      - Series:   9 {ep2d_diff, mb3, ref                       } [  |    |    45] 
      - Series:  10 {ep2d_diff, mb3, blip inv                  } [  |    |    45] 
      - Series:  11 {mprage, HR64                              } [  |    |   240] 
      - Series:  12 {localizer (expectb-b, blank)              } [  |    |    11] 
      - Series:  13 {AAHead_Scout_64ch-head-coil               } [  |    |   128] 
      - Series:  14 {AAHead_Scout_64ch-head-coil_MPR_sag       } [  |    |     5] 
      - Series:  15 {AAHead_Scout_64ch-head-coil_MPR_cor       } [  |    |     3] 
      - Series:  16 {AAHead_Scout_64ch-head-coil_MPR_tra       } [  |    |     3] 
      - Series:  17 {gre_field_map, 2mm, filter M              } [  |    |    90] 
      - Series:  18 {gre_field_map, 2mm, filter M              } [  |    |    45] 
      - Series:  19 {ep2d_bold, mb3, part 1                    } [  |    |  1061] 
      - Series:  20 {ep2d_bold, mb3, part 2                    } [  |    |   887] 
      - Series:  21 {ep2d_bold, mb3, part 3                    } [  |    |  1061] 
      - Series:  22 {ep2d_diff, mb3, ref                       } [  |    |    45] 
      - Series:  23 {ep2d_diff, mb3, blip inv                  } [  |    |    45] 
=====
The latest recorded HR data:
Series:  11 {mprage, HR64                              } [  |    |   240] 
DicomDownload:
Calling system's COPY function to dump the data...16:05:04
source:/XXX/XXX/prisma/images/XXXX
destination:mrt/
COPY finished successully 16:05:04
ConvertDicom:
Found 240 files...
Dicom conversion s#1... (16:05:10)
Will call spm_jobman...
Running SPM jobman 1...
------------------------------------------------------------------------
Running job #1
------------------------------------------------------------------------
Running 'DICOM Import'
   Changing directory to: mrt/
   Changing back to directory: /tmp/fancycarp
Done    'DICOM Import'
Done
Finished... (16:05:15)
Deleting DICOM images in (16:05:15)
mrt/
Finished... (16:05:15)

6/ Downloading Functional Data

Similarly, get_epi method downloads the EPIs based on ProjectObject properties.

>> s.get_epi
Making a dicom query, sometimes this might take long (so be patient)...(16:05:14)
This is what I found for you:
You told me to download the following series: 8,19,20,21,6,7,17,18,
Double check if everything is fine.
Making a dicom query, sometimes this might take long (so be patient)...(16:05:15)
Will now dump series (16:05:17)
DicomDownload:
Calling system's COPY function to dump the data...16:05:17
source:/XXX/XXX
destination:/tmp/mynextsciencepaper/data//sub001/run001/mrt/

7/ Preprocessing

s.preprocess_pipeline takes care of the preprocessing steps. It tries to make a fieldmap correction if required files are present. And continues with a/ Surface Segmentation of the anatomical data, b/ Normalization to MNI space, c/ Realignment time and participants, d/ Gray to white matter segmentation with newSegment, e/ Strips away skull voxels

 function preprocess_pipeline(self,runs)
            %meta method to run all the required steps for hr
            %preprocessing. RUNS specifies the functional runs, make it a
            %vector if needed. RUNS will be used with Re_Coreg.
            if nargin > 1
                try
     fprintf('Will attempt field corrections\n');
     self.ComputeVDM;
                 self.ApplyVDM;
     self.epi_prefix = 'u_';%if we came that far, change the default EPI prefix.
    catch
     fprintf('Failed... Will work on non-field corrected EPIs\n');
    end
                %
       self.SegmentSurface_HR;%cat12 segmentation
             self.SkullStrip;%removes non-neural voxels
             self.MNI2Native;%brings the atlas (if present) to native space
             self.Re_Coreg(runs);%realignment and coregistration
             self.Segment_meanEPI;%segments mean EPI with new segment
    self.SkullStrip_meanEPI;%creates a native mask
      else
    fprintf('One input argument is required!\n');
      end
        end

8/ Analysis of functional data. So far, Fancycarp covered all the basic preprocessing steps, which should be fairly common to anybody, independent of their projects. However, first-level analyses in fMRI has certainly project specific flavors, that cannot be fully automated with additional data. To this end, the CreateFolderHierarchy method has created the design folder in every subject and run folder, which is supposed to contain parameters necessary for the building of a design matrix required for the first-level analysis. One may create as many design models indexed with an integer, and use this number as argument to first-level analysis method defined in the SubjectObject.

>> pwd
ans =
/mnt/data/project_fearamy/data/sub005/run001/design
>> ls -l
total 136
drwxr-xr-x 2 onat onat 4096 Feb 24  2017 model01
drwxr-xr-x 2 onat onat 4096 Feb 24  2017 model02
drwxrwxr-x 2 onat onat 4096 Feb 13 11:29 model03
drwxr-xr-x 2 onat onat 4096 Feb 24  2017 model04
....

In each of these folders there should be a data.mat file that contains SPM compatible timeing information for different stimuli and conditions. For example the model below contains 10 condition as a Matlab structure which could be directly fed to the SPM gui.

>> cond
cond = 
1x10 struct array with fields:
    name
    onset
    duration
    tmod
    pmod

Calling s.analysis_firstlevel(1,1) should fit a first-level GLM using the model 1 to run 1.

Guide for reproducing our BioRxiv Manuscript.

In this post, I would like to provide guidelines on how to reproduce results and figures presented in our recent manuscript in BioRxiv. I will show how to download the raw data, run the analysis routines,  and produce figures that are shown in the paper. The text below is also present in the READ.me section of the associated Matlab toolbox that is published as an open-source in GitHub.

About

With this repository you can reproduce all the analysis and figures presented in our paper:

Fear Generalization as Threat Prediction: Adaptive Changes in Facial Exploration Strategies revealed by Fixation-Pattern Similarity Analysis. Bioarxiv

In this paper, we investigate how humans explore faces and focus how these exploration patterns change when faces are associated with an aversive outcome.

Initial Setup

You can download the data and the associated repositories, including this one from the Open Science Framework here. Add these repositories to your Matlab path.

addpath('/home/onat/Documents/Code/Matlab/FPSA_FearGen/');
addpath('/home/onat/Documents/Code/Matlab/globalfunctions//');
addpath('/home/onat/Desktop/fancycarp');

Examples

Get the list of participants.

The data you downloaded from the OSF contains all recorded participants. However, few had to be excluded for the main analysis. It is important to include all participants as it provide the possibility to reproduce the selection criteria used in the report or test the results with another selection criteria. 'get_subjects' action returns all the selected participants, totalling to 74.

>> FPSA_FearGen('get_subjects')
ans =
  Columns 1 through 45
     1     2     3     4     5     6     7     8     9    10    11    12    14    15    16    17    18    19    20    21    22    23    24    25    26    27    28    29    30    31    32    33    35    36    37    39    40    41    42    43    44    45    46    47    48
  Columns 46 through 74
    49    50    51    53    54    56    57    59    60    61    62    63    64    65    66    67    68    69    70    71    72    73    74    77    78    79    80    81    82
>> 

Sanity check 1.

In eye tracking it is usual that few trials are excluded due to blinks or bad calibration. The following code plots the number of trials per condition across the participant pool. The figure shows that almost all the particpants have their trials recorded as expected. This sanity check also shows that few participants have less trials than others.

FPSA_FearGen('get_trialcount',4)

4 above refers to the generalization phase. (2: baseline; 3: conditioning; 4: test phase).

Get the data to workspace using the get_fixmat action.

Fixmat is a compact way of storing large amounts of eye movement recordings in the form of fixation points.

>> fixmat = FPSA_FearGen('get_fixmat')
ans = 
  Fixmat with properties:

     subject: [1×118188 uint32]
       phase: [1×118188 int32]
       start: [1×118188 int32]
        stop: [1×118188 int32]
           x: [1×118188 single]
           y: [1×118188 single]
         eye: [1×118188 uint8]
    deltacsp: [1×118188 int32]
        file: [1×118188 int32]
     oddball: [1×118188 int32]
     trialid: [1×118188 int32]
         ucs: [1×118188 int32]
         fix: [1×118188 int32]
       chain: [1×118188 double]
       isref: [1×118188 double]

It stores every fixation's attributes in the form of separate vectors. For example, information such as x and y coordinates, participant's index, the image, the condition are stored in the Fixmat.

However, Fixmat variable above is not a simple Matlab structure though, but rather an instance of a Fixmat object, as coded in the FancyCarp Toolbox. The benefit of having a Fixmat object is that there are useful methods built-in in the Fixmat object, such as for example visualizing fixation density maps as heatmaps.

For example, the following code can be used to plot a fixation density map (FDM) based on all subjects (stored in S) during the generalization phase:

S = FPSA_FearGen('get_subjects');
fixmat = FPSA_FearGen('get_fixmat');
v{1} = {'subject' S 'phase' 4};
fixmat.getmaps(v{:});
fixmat.plot

The Fixmat.getmaps method creates an FDM based on the cell array argument v. In the example above,v is used to select all fixations that belong to both phase 4 and subjects S. The method Fixmat.plot plots the computed FDM. In order to create a separate FDM for different conditions or participants we would create a separate cell array for each of the required filtering conditions.

v=[];
c=0;
for ns = S([11 3 8 23]);
  c=c+1;
  v{c} = {'subject' ns 'deltacsp' 0 'phase' 4};
end
fixmat.getmaps(v{:});
fixmat.plot

This will plot a separate FDM for the 4 different participants.

Note how different participants scan faces differently. This is the basis for recent reports that analyzed the scan-path idiosyncrasy during viewing of faces. And the major reason for us to come up with the FPSA as a methodology to investigate how eye movement strategies change with aversive learning during viewing of faces.

Get FDM for the Fixation-Pattern Similarity Analysis

FPSA analysis is conducted on single-participants at a time. The get_fixmap action can be used to gather the required FDMs. It is basically a wrapper around the Fixmat.getmaps method.

>> maps = FPSA_FearGen('get_fixmap',fixmat,{'subject' 2});;
>> size(maps)
ans =
      250000          16

As you can see in maps each FDM is represented as a column vector. As we have 2 phases and 8 different conditions (faces), this amounts to 16 different vectors. This representation is appropriate for the similarity analysis as it can be readily used as an argument to Matlab's pdist

Fixation-Pattern Similarity Analysis

The following command will run a correlation-based similarity analysis on the FDMs of single participants, before (phase = 2) and after (phase =4) aversive learning. {'fix',1:100} indicates that the analysis will include 1st to 100th fixations, that is all the fixations.

1:3 ensures that similarity analysis is ran separately for the 3 runs of the generalization phase (phase = 4). This is important because the baseline phase (phase = 2) entails only one single run. In order to have a valid comparison of the similarity values across the baseline and test phases, the FDMs should contain on average same number of fixations. Computing separate FPSA matrices for each run ensures this.

sim = FPSA_FearGen('get_fpsa_fair',{'fix',1:100},1:3);
>> sim
sim = 
  struct with fields:

    correlation: [74×120 double]

The resulting matrix sim.correlation contains in each row the similarity data from a given participant. The pair-wise similarity values between the 16 FDMs are stored across the columns in a non-redundant format. You can visualize the resulting 16x16 matrix after putting the similarity values back to their positive-definite matrix format with the squareform function. 'plot_fpsa' action is used to do this.

The first quadrant of this matrix shows the similarity relationships between the 8 fixation-density maps recorded during the baseline period, before aversive learning has taken place. The second quadrant shows the same after learning has taken place.

Multi-dimensional Similarity Analysis

Another extremely useful way of representing the similarity relationships between the fixation maps consists of using the multi-dimensional scaling method. Below, I ran an MDS analysis using two dimensions.

fixmat = FPSA_FearGen('get_mdscale',squareform(mean(sim.correlation)),2)

This places a set of dots in such a way that the distances between each dot corresponds to the dissimilarity between the FDMs. Therefore, it is nice visual tool that summarizes the complex disssimilarity matrice, which is sometimes hard to digest. 

With the following piece of code, the upper part of the figure 03 of the manuscript can be produced summing up bundling previous figures.

FPSA_FearGen('figure_03A',FPSA_FearGen('get_fpsa_fair',{'fix' 1:100},1:3))

Analysis of dissimilarity matrices

One major aim in our paper was to understand how the similarity relationships between the FDMs changed with aversive learning. We take a linear modelling approach:

Y = bX

where Y is a non-redundant dissimilarity matrix, X is one of the 3 models we are testing, b are the coefficients that are to be fitted.

>> FPSA_FearGen('FPSA_get_table',{})
ans = 
      FPSA_B        FPSA_G         circle        specific      unspecific     Gaussian    subject    phase
    __________    ___________    ___________    ___________    ___________    ________    _______    _____
      -0.52828      -0.014812        0.70711     4.3298e-17        0.70711    0.48254     1          1    
     -0.064526        -0.3752    -1.0147e-17           -0.5            0.5    0.32924     1          1    
      -0.21919        -0.3426       -0.70711       -0.70711     8.6596e-17    0.25385     1          1    
      -0.28268       -0.20612             -1           -0.5           -0.5    0.32924     1          1    
      -0.35528      -0.052057       -0.70711    -1.2989e-16       -0.70711    0.48254     1          1    
       0.10652     -0.0068816    -1.7938e-16            0.5           -0.5    0.50102     1          1    
       0.27702        0.12585        0.70711        0.70711    -1.7319e-16    0.49959     1          1    
      -0.05106       -0.13439        0.70711    -4.3298e-17        0.70711    0.56373     1          1    
      -0.10373      -0.067047     6.1232e-17    -6.1232e-17     1.2246e-16      0.523     1          1    
       -0.1175       -0.39994       -0.70711    -4.3298e-17       -0.70711    0.56373     1          1  
       ...

The first two columns are the observed dissimilarity values for the first subject in the first phase (see last columns). Columns with names circle, specific, unspecific and Gaussian are the modelled dissimilarity values that will be fitted for each participant separately.

In Matlab's statistical toolbox one can fit a linear model (pretty much the same way as in R) as follows:

fitlm(T,'FPSA_B ~ 1 + circle');

To to this for all the participants and models separately, we can run

T = FPSA_FearGen('FPSA_get_table',{});
[A,C] = FPSA_FearGen('FPSA_model_singlesubject',T)  
>> A.model_01
ans = 
  struct with fields:

    w1: [74×2 double]
>> A.model_01
ans = 
  struct with fields:

    w1: [74×2 double]
>> A.model_02
ans = 
  struct with fields:

    w1: [74×2 double]
    w2: [74×2 double]

The results of this analysis can be plotted with the following code, reproducing bottom panels of the figure 03.

FPSA_FearGen('figure_03C',{'fix' 1:100})