Dear Michael,

Sorry about the delay in replying.
I have double-checked with the code and can now give you
definitive answers to your questions.

Possum simulates each tissue separately and uses the
fractional occupancy value at each voxel to weight the
signal contribution from that voxel.  So there is no hard
classification involved - each tissue's signal is weighted by 
the partial volume (or occupancy) value, which is a floating
point number between 0 and 1.

What you do about the T2* changes is purely up to you,
but the way that possum models things is to apply the
specified T2* changes equally to _all_ tissues present
in the voxel.  You can arrange for the change to be zero
in voxels that do not contain GM, or you can include it.
This choice is up to you and what you want the simulation
to reflect.  

If you wanted to only have T2* changes associated
with GM then you can either just zero the T2* changes
in any voxel that has no GM and then leave the ones with a 
fraction of GM (i) at their full value or (ii) scale the T2* changes
with the GM fraction.  Alternatively, you could split up the
input object into three separate tissues (rather than a
4D object containing all three) and then do three separate
simulations with the ability to be completely free to do whatever
you want with each tissue totally separately.  I don't know if it
will make a lot of difference, or even what the "correct" way is
to simulate these changes in voxels where there is a mixture
of GM and either WM or CSF.  If you have a strong theory
on this then I hope that one of the above options will give
you what you want.

All the best,

On 2 Feb 2012, at 22:11, Michael Hallquist wrote:

> Hi Mark,
> I'm not sure I understand. I see that the input for possum is
> segmented by tissue type, but the data are fractional and appear to
> represent the proportion (or perhaps probability) of a tissue type at
> that location. How does POSSUM handle the situation, for example,
> where a voxel has a value in the GM map of 0.4 and the WM map is 0.6?
> I notice in the POSSUM logs that the core voxelwise loop iterates from
> tissue type 0-2, corresponding to GM, WM, and CSF. Does POSSUM do a
> "hard classification" of the 0-1 values in the tissue segmentation
> maps?
> More generally, if I am trying to simulate activations within a
> spherical ROI where the sphere includes 515 voxels (10mm diameter),
> then within the sphere there are voxels that are predominantly WM, GM,
> or CSF. Should I include activation time courses (even if these are
> relatively flat given that T2* changes are predominant in GM tissue in
> the BOLD contrast) in these voxels, or should these be left at 0?
> When you say everything else will be taken care of by the simulation,
> do you mean that if a voxel contains 0.3 WM, for example, along with
> 0.7 GM, and I provide a corresponding activation time course (computed
> as change against the static T2* value for GM), then the activation
> time course will only influence that part of signal intensity
> attributable to GM, whereas the 0.3 WM will influence signal intensity
> only through its static T2* value? Lastly, what do you mean about
> running the simulations for each tissue type separately? Is this
> something I would do manually or are you referring to the GM, WM, CSF
> input to POSSUM and superordinate loop over tissue type?
> Sorry for the blitz of questions. I want to be careful to get this
> right since the underlying rationale of the simulation depends on
> accurate extraction of the activation time courses.
> Thanks again,
> Michael
> On Thu, Feb 2, 2012 at 3:05 PM, Mark Jenkinson <[log in to unmask]> wrote:
>> Hi Michael,
>> It is actually simpler than this.
>> There is a separate input for GM, WM and CSF (look at the brain.nii.gz
>> image in $FSLDIR/data/possum/).  So what we recommend is that you
>> run the simulations for each tissue type separately.  Therefore you only
>> need to specify a changing T2* for the GM tissue type.  Everything else
>> will be taken care of by the simulations.
>> As for masking/ROI, I think you're strategy sounds fine, although you can
>> compare some small scale simulations with and without this masking/ROI
>> in order to double check.
>> All the best,
>>        Mark
>> On 2 Feb 2012, at 17:26, Michael Hallquist wrote:
>>> Dear Mark,
>>> Many thanks for your detailed response! I now have a good handle on
>>> this issue and have been able to setup the 4d activation files
>>> successfully.
>>> I had a couple of quick follow-up questions. I am placing activation
>>> time courses within spherical ROIs for the simulation and the spheres
>>> often include somewhat different tissue types. I've used FAST with a
>>> single subject's brain to get the GM/WM/CSF segmentation values, which
>>> sum to 1. Is a reasonable approach for scaling a time series in a
>>> voxel that has, for example, 0.61 GM and 0.39 WM, to use a linear
>>> combination based on the static T2* values? Using your TE=30ms 10%
>>> signal change example, would the proper scaling be:
>>> [0.61*[(30/(30/51-ln(1.1))) - 51]+ 0.39*[(30/(30/44-ln(1.1))) - 44]]/1000
>>> (sorry that's so ugly to write by email)
>>> In terms of the masking of brain regions to speed up the simulation, I
>>> am interested in the effects of motion, so it sounds problematic not
>>> to simulate the whole volume. That said, if I were to start with a
>>> template that has activations in just two 10mm spherical ROIs with a
>>> maximum motion displacement of 1mm, would I be able to obtain accurate
>>> estimates by masking the segmented input file using the ROIs, but
>>> dilating the mask 5 or 6mm to ensure that even with motion, the ROI
>>> would be properly sampled?
>>> Thanks,
>>> Michael
>>>> You can limit the voxels that POSSUM calculates by changing
>>>> the input brain object.  It will only calculate things for object
>>>> voxels present in this image.  So feel free to take the default
>>>> brain and mask out whatever you like.  This will certainly speed
>>>> things up, although it won't give you the right solutions in the
>>>> case of motions, as you need the other voxels present, but I
>>>> don't know if that is of interest to you or not.
>>> On Wed, Jan 25, 2012 at 2:36 PM, Mark Jenkinson <[log in to unmask]> wrote:
>>>> Dear Michael,
>>>> You've pretty much got it.
>>>> The list of points 1-5 are all correct except for:
>>>>  - units are in seconds, not milliseconds
>>>>  - T2* changes are all relative to the baseline, static value
>>>>        and not the prior volume (as 1 seems to imply, but
>>>>        conflicts with what you say in 2, which is correct)
>>>> So you need to scale the values in your 4D file to represent
>>>> changes in T2* in units of seconds (e.g. 0.001 would represent
>>>> 1 millisecond).
>>>> S0 is an intensity value not a T2* value, so it is in arbitrary units
>>>> and not in seconds or milliseconds.  The S = S0*exp(-TE/T2*)
>>>> equation is about _intensity_ values.  So if you had a TE of 30ms
>>>> and a T2* of 51ms then S=S0*exp(-30/51) = 0.55531*S0
>>>> whereas if T2* increased by 5.1ms (an increase is what you'd
>>>> expect with more activation, as there's less deoxyHb and less
>>>> field distortion, so longer T2*) then you'd have an intensity of
>>>> S=S0*(-30/56.5) = 0.58803*S0.  This would represent a
>>>> signal increase of 0.58803/0.55531 = 1.0589 (or about 6%).
>>>> Note that this signal increase is more simply calculated as
>>>> exp(-TE/T2*B +TE/T2*1) where T2*A and T2*B are the two relevant
>>>> T2* values (so exp(-30/56.5 +30/51)=1.0589).
>>>> Therefore, to achieve a certain value of percent signal change
>>>> (which is normally what people think about) then you need to
>>>> solve the equation:  exp(-TE/T2*B + TE/T2*A)=desired factor
>>>> setting T2*A=51ms and TE and desired factor as appropriate.
>>>> For example, a 10% signal change, with a 30ms TE, means:
>>>>  exp(-30/T2* + 30/51)=1.1  => -30/T2* + 30/51=ln(1.1)
>>>>  =>  T2* = 30/(30/51-ln(1.1)) = 60.8 and so the change in T2*
>>>> that would need to be encoded in the 4D activation file would
>>>> be 60.8 - 51 = 9.8ms = 0.0098 seconds.
>>>> I hope this helps explain how this part works.
>>>> You do need to specify the --activt4D input.
>>>> You can limit the voxels that POSSUM calculates by changing
>>>> the input brain object.  It will only calculate things for object
>>>> voxels present in this image.  So feel free to take the default
>>>> brain and mask out whatever you like.  This will certainly speed
>>>> things up, although it won't give you the right solutions in the
>>>> case of motions, as you need the other voxels present, but I
>>>> don't know if that is of interest to you or not.
>>>> All the best,
>>>>        Mark
>>>> On 24 Jan 2012, at 19:37, Michael Hallquist wrote:
>>>>> Dear FSL gurus (especially Mark and Ivana):
>>>>> I have a question about how to setup a 4d activation timecourse file
>>>>> using simulated time series not generated from empirical fMRI data. I
>>>>> realize that there have been a number of questions on the listserv
>>>>> about how to specify these files and I've read through these posts in
>>>>> detail. Yet I'm getting stuck on how to properly scale time series
>>>>> that are in arbitrary units and where there is no task baseline per se
>>>>> (resting state data).
>>>>> Here's the summary of what I've understood (and please let me know if
>>>>> I've misunderstood!):
>>>>> 1) The 4d activation file provides voxelwise timecourse specifications
>>>>> of T2* changes (i.e., units are T2star dt [delta time], relative to
>>>>> prior volume).
>>>>> 2) The time course units are milliseconds (representing relaxation
>>>>> time) and are relative to the baseline T2* relaxation value specified
>>>>> in the MR parameters file. In the case of the stock MRpar_3T file, the
>>>>> T2* relaxation time for  gray matter is 51ms.
>>>>> 3) A value of 0 at any point along a voxel timecourse in the 4d
>>>>> activation file means that the intensity will be based on the static
>>>>> T2* value alone (51ms per above).
>>>>> 4) T2* values are related to simulated signal intensity by the
>>>>> following equation: S = S0 * exp(-TE/T2*), where S0 is the baseline
>>>>> intensity, TE is the echo time, and T2* is the absolute T2* value.
>>>>> 5) To generate signal intensity at a given voxel and timepoint, POSSUM
>>>>> sums the T2* value of the static tissue type (51ms above) plus the
>>>>> voxelwise, time-specific value in the 4d activation file (e.g., 5.1).
>>>>> Now to my questions. I have simulated correlated time series in
>>>>> arbitrary units that are scaled roughly [-2, 2]. How do I scale these
>>>>> time series into the appropriate T2* units for the --activ4D input
>>>>> given the equation above? In my case, is S0 51ms or do I need some
>>>>> baseline intensity in arbitrary scanner units?
>>>>> Also, If the static T2* time for gray matter is 51ms and I provide a
>>>>> value of 5.1 at the given point in the timecourse, does that
>>>>> correspond to a 10% signal change?
>>>>> Do I always need to specify --activt4D to provide the times for the 4d
>>>>> activation file, or will these be computed based on the TR encoded in
>>>>> the NIfTI of the activation file?
>>>>> Lastly and most peripherally, I realize that the 4d activation results
>>>>> in computationally and RAM-intensive simulations, which is fine. Is
>>>>> there a way to reduce the RAM or computation burden using a mask file
>>>>> (simulation within a sparse matrix of sorts)?
>>>>> Thanks very much for your help,
>>>>> Michael