diff --git a/Examples/Scripts/IcaFixProcessingBatch.sh b/Examples/Scripts/IcaFixProcessingBatch.sh
index c955c6785..d081eaeba 100755
--- a/Examples/Scripts/IcaFixProcessingBatch.sh
+++ b/Examples/Scripts/IcaFixProcessingBatch.sh
@@ -153,6 +153,9 @@ main() {
 	# set whether or not to regress motion parameters (24 regressors)
 	# out of the data as part of FIX (TRUE or FALSE)
 	domot=FALSE
+
+	# set the ICA Method (MELODIC or ICASSO)
+	ICAmethod=MELODIC
 	
 	# set the training data used in multi-run fix mode
 	MRTrainingData=HCP_Style_Single_Multirun_Dedrift.RData
@@ -233,7 +236,7 @@ main() {
 
 				echo "  InputFile: ${InputFile}"
 
-				cmd=("${queuing_command[@]}" "${HCPPIPEDIR}/ICAFIX/hcp_fix_multi_run" --fmri-names="${InputFile}" --high-pass=${bandpass} --concat-fmri-name="${ConcatFileName}" --motion-regression=${domot} --training-file="${MRTrainingData}" --fix-threshold=${FixThreshold} --delete-intermediates="${DeleteIntermediates}" --config="$config" --processing-mode="$processingmode")
+				cmd=("${queuing_command[@]}" "${HCPPIPEDIR}/ICAFIX/hcp_fix_multi_run" --fmri-names="${InputFile}" --high-pass=${bandpass} --concat-fmri-name="${ConcatFileName}" --motion-regression=${domot} --training-file="${MRTrainingData}" --fix-threshold=${FixThreshold} --delete-intermediates="${DeleteIntermediates}" --config="$config" --processing-mode="$processingmode" --ica-method="$ICAmethod")
 				echo "About to run: ${cmd[*]}"
 				"${cmd[@]}"
 			done
diff --git a/ICAFIX/hcp_fix_multi_run b/ICAFIX/hcp_fix_multi_run
index 9792a7ca7..dc7201b5d 100755
--- a/ICAFIX/hcp_fix_multi_run
+++ b/ICAFIX/hcp_fix_multi_run
@@ -28,6 +28,10 @@
 #   1) Improvements to script readability and consistency
 #   2) Addition of a polynomial detrend option
 #   3) Addition of FIX threshold as optional parameter
+#
+#   Changes by Burke Rosen
+#   1) Optionally use icasso instead of melodic for IC estimation and use melodic for mixture modeling only
+#      (Followed by FIX as usual)
 
 #############################################################
 set -eu
@@ -126,6 +130,13 @@ opts_AddOptional '--processing-mode' 'ProcessingMode' '"HCPStyleData" (default)
 
 opts_AddOptional '--concatenate-only' 'ConcatOnlyStr' 'TRUE or FALSE' "When set to TRUE the script stops after the concatination step. e.g., for use in experimental alternative multi-run denoising" 'FALSE'
 
+opts_AddOptional '--ica-method' 'ICAmethod' 'MELODIC or ICASSO' "Use single-pass MELODIC (default) or multi-pass ICASSO consensus method for ICA" 'MELODIC'
+
+opts_AddOptional '--icasso-vis' 'icasso_vis' 'basic or off' "Save icasso figures (default) or do not" 'basic'
+
+opts_AddOptional '--icasso-nICA' 'icasso_nICA' '100' "Number of ICA repetitions per icasso repetition, @ delimited string (default = 100)" 'basic'
+
+
 opts_ParseArguments "$@"
 
 if ((pipedirguessed))
@@ -287,20 +298,20 @@ demeanMovementRegressors() {
 	AllOut=""
 	c=1
 	while (( c <= nCols )) ; do
-		ColIn=`cat ${In} | sed 's/  */ /g' | sed 's/^ //g' | cut -d " " -f ${c}`
+		ColIn=$(cat ${In} | sed 's/  */ /g' | sed 's/^ //g' | cut -d " " -f ${c})
 		bcstring=$(echo "$ColIn" | tr '\n' '+' | sed 's/\+*$//g')
 		valsum=$(echo "$bcstring" | bc -l)
 		valmean=$(echo "$valsum / $nRows" | bc -l)
 		ColOut=""
 		r=1
 		while (( r <= nRows )) ; do
-			val=`echo "${ColIn}" | head -${r} | tail -1`
-			newval=`echo "${val} - ${valmean}" | bc -l`
-			ColOut=`echo ${ColOut} $(printf "%10.6f" $newval)`
+			val=$(echo "${ColIn}" | head -${r} | tail -1)
+			newval=$(echo "${val} - ${valmean}" | bc -l)
+			ColOut=$(echo ${ColOut} $(printf "%10.6f" $newval))
 			r=$((r+1))
 		done
-		ColOut=`echo ${ColOut} | tr ' ' '\n'`
-		AllOut=`paste <(echo "${AllOut}") <(echo "${ColOut}")`
+		ColOut=$(echo ${ColOut} | tr ' ' '\n')
+		AllOut=$(paste <(echo "${AllOut}") <(echo "${ColOut}"))
 		c=$((c+1))
 	done
 	echo "${AllOut}" > ${Out}
@@ -352,6 +363,17 @@ then
 	log_Err_Abort "--fallback-threshold is not allowed when --vol-wisharts=1 is used"
 fi
 
+case "$ICAmethod" in
+	(MELODIC|ICASSO)
+		;;
+	(*)
+		log_Err_Abort "--ica-method=${ICAmethod} must be MELODIC or ICASSO"
+		;;
+esac
+if [[ ${FSL_FIX_MATLAB_MODE} == 0 ]] && [[ ${ICAmethod} == ICASSO ]];then
+	log_Err_Abort "--ica-method=ICASSO is not yet implemented for compiled matlab."
+fi
+
 ## ---------------------------------------------------------------------------
 ## Report parameters
 ## ---------------------------------------------------------------------------
@@ -366,7 +388,7 @@ then
 	log_Warn "--icadim-mode='fewtimepoints' is being used with multiple wisharts, multiple wishart fitting is not expected to work well when the data has few timepoints"
 fi
 
-DIR=`pwd`
+DIR=$(pwd)
 log_Msg "PWD: $DIR"
 
 # For INTERPRETED MODES, make sure that matlab/octave has access to the functions it needs.
@@ -420,11 +442,11 @@ for fmri in "${fmris[@]}" ; do
 	MovementNIFTIMergeARRAY+=("${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf")
 	MovementNIFTIhpMergeARRAY+=("${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp")
 	
-	cd `dirname $fmri`
+	cd $(dirname $fmri)
 	log_Debug_Msg "pwd: "$(pwd)
-	fmri=`basename $fmri`  # After this, $fmri no longer includes the leading directory components
-	fmri=`$FSLDIR/bin/imglob $fmri`  # After this, $fmri will no longer have an extension (if there was one initially)
-	if [ `$FSLDIR/bin/imtest $fmri` != 1 ]; then
+	fmri=$(basename $fmri)  # After this, $fmri no longer includes the leading directory components
+	fmri=$($FSLDIR/bin/imglob $fmri)  # After this, $fmri will no longer have an extension (if there was one initially)
+	if [ $($FSLDIR/bin/imtest $fmri) != 1 ]; then
 		log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}"
 	fi
 
@@ -451,7 +473,7 @@ for fmri in "${fmris[@]}" ; do
 	
 	## "1st pass" VN on the individual runs; high-pass gets done here as well
 	#fslval puts a space on the end of its number, remove it...
-	tr=`$FSLDIR/bin/fslval $fmri pixdim4 | tr -d ' '`  #No checking currently that TR is same across runs
+	tr=$($FSLDIR/bin/fslval $fmri pixdim4 | tr -d ' ')  #No checking currently that TR is same across runs
 	log_Msg "tr: ${tr}"  
 	log_Msg "processing FMRI file ${fmri} with highpass ${hp}"
 
@@ -511,7 +533,6 @@ for fmri in "${fmris[@]}" ; do
 
 done  ### END LOOP (for fmri in $fmris; do)
 
-
 #note: this test a string comparison, not a math comparison
 if [[ "$RewishartThresh" != "0" ]]
 then
@@ -540,10 +561,10 @@ then
 			#rerun with 1 wishart
 			fmri="${fmriarray[$i]}"
 			#fmri is modified by basename and some other stuff, tr is derived from file...
-			cd `dirname $fmri`
+			cd $(dirname $fmri)
 			log_Debug_Msg "pwd: "$(pwd)
-			fmri=`basename $fmri`  # After this, $fmri no longer includes the leading directory components
-			fmri=`$FSLDIR/bin/imglob $fmri`  # After this, $fmri will no longer have an extension (if there was one initially)
+			fmri=$(basename $fmri)  # After this, $fmri no longer includes the leading directory components
+			fmri=$($FSLDIR/bin/imglob $fmri)  # After this, $fmri will no longer have an extension (if there was one initially)
 
 			#log the fact that we are rerunning this one
 			log_Msg "rerunning functionhighpassandvariancenormalize on outlier run $fmri"
@@ -551,11 +572,11 @@ then
 			#TSC: assume _mean and _demean haven't been overwritten or deleted?
 			
 			## "1st pass" VN on the individual runs; high-pass gets done here as well
-			tr=`$FSLDIR/bin/fslval $fmri pixdim4 | tr -d ' '`  #No checking currently that TR is same across runs
+			tr=$($FSLDIR/bin/fslval $fmri pixdim4 | tr -d ' ')  #No checking currently that TR is same across runs
 			log_Msg "tr: ${tr}"  
 			log_Msg "processing FMRI file ${fmri} with highpass ${hp}"
 
-			if [ `$FSLDIR/bin/imtest $fmri` != 1 ]; then
+			if [ $($FSLDIR/bin/imtest $fmri) != 1 ]; then
 				log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}"
 			fi
 
@@ -608,7 +629,7 @@ fi
 ## ---------------------------------------------------------------------------
 
 #Make Concatenated Folder
-ConcatFolder=`dirname ${ConcatName}`
+ConcatFolder=$(dirname ${ConcatName})
 log_Msg "Making concatenated folder: ${ConcatFolder}"
 if [ ! -e ${ConcatFolder} ] ; then  
 	mkdir ${ConcatFolder} 
@@ -689,13 +710,13 @@ if ((ConcatOnly)); then
 fi
 
 ## ---------------------------------------------------------------------------
-## Prepare for melodic on concatenated file
+## Prepare for melodic/icasso on concatenated file
 ## ---------------------------------------------------------------------------
 
 # concatfmri is intensity normalized and has no overall mean.
 # It is also not variance normalized (the mean VN map *across runs* was "restored" above).
 # But is detrended and within run demeaned and variance normalized for both volume and CIFTI
-concatfmri=`basename ${ConcatNameNoExt}`  # Directory path is now removed
+concatfmri=$(basename ${ConcatNameNoExt})  # Directory path is now removed
 concatfmrihp=${concatfmri}_hp${hp}
 
 #MPH: Need to be in original directory here, for the ${MovementNIFTI*ARRAY[@]} variables to work
@@ -740,79 +761,93 @@ case ${FSL_FIX_MATLAB_MODE} in
 		;;
 esac
 
-Dim=`cat ${concatfmrihp}_dims.txt`
+Dim=$(cat ${concatfmrihp}_dims.txt)
 
 log_Msg "Dim: ${Dim}"
 
-## ---------------------------------------------------------------------------
-## Run melodic on concatenated file
-## ---------------------------------------------------------------------------
+if [[ $ICAmethod == MELODIC ]]; then
+	## ---------------------------------------------------------------------------
+	## Run melodic on concatenated file
+	## ---------------------------------------------------------------------------
 
-log_Debug_Msg "About to run melodic: Contents of ${concatfmrihp}.ica follow"
-if [ ! -z "${log_debugOn}" ] ; then
-	ls -lRa ${concatfmrihp}.ica
-fi
+	log_Debug_Msg "About to run melodic: Contents of ${concatfmrihp}.ica follow"
+	if [ ! -z "${log_debugOn}" ] ; then
+		ls -lRa ${concatfmrihp}.ica
+	fi
 
-#grab melodic from $PATH by default, don't hardcode it with respect to $FSLDIR
-#we need to do "if which ..." because the script currently uses ERR trap
-if which melodic &> /dev/null
-then
-	MELODIC=$(which melodic 2> /dev/null)
-else
-	#if it isn't even in $PATH, fall back on FSLDIR
-	MELODIC="${FSLDIR}/bin/melodic"
-fi
+	#grab melodic from $PATH by default, don't hardcode it with respect to $FSLDIR
+	#we need to do "if which ..." because the script currently uses ERR trap
+	if which melodic &> /dev/null
+	then
+		MELODIC=$(which melodic 2> /dev/null)
+	else
+		#if it isn't even in $PATH, fall back on FSLDIR
+		MELODIC="${FSLDIR}/bin/melodic"
+	fi
 
-log_Msg "Running MELODIC located at: $MELODIC"
-log_Debug_Msg "Beginning of melodic version log, help, and checksum"
-if [ ! -z "${log_debugOn}" ] ; then
-	log_Debug_Msg "$MELODIC --version"
-	$MELODIC --version
-	log_Debug_Msg "$MELODIC --help"
-	$MELODIC --help
-	log_Debug_Msg "md5sum $MELODIC"
-	md5sum $MELODIC
-fi
-log_Debug_Msg "End of melodic version log, help, and checksum"
-
-# Use the Wishart-based dim estimate instead of the melodic-based estimate (--dim=${Dim})
-# Also turn off melodic VN and feed in the input that has already been VN'ed instead (--vn flag, which turns OFF melodic's vn)
-# Added "-m ${concatfmri}_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018
-melodic_cmd=("${MELODIC}" -i "${concatfmrihp}_vnts" -o "${concatfmrihp}.ica/filtered_func_data.ica" --nobet --report --Oall --tr="${tr}" --vn --dim="${Dim}" -m "${concatfmri}_brain_mask")
-if [ ! -z "${log_debugOn}" ] ; then
-	melodic_cmd+=(--verbose --debug)
-fi
+	log_Msg "Running MELODIC located at: $MELODIC"
+	log_Debug_Msg "Beginning of melodic version log, help, and checksum"
+	if [ ! -z "${log_debugOn}" ] ; then
+		log_Debug_Msg "$MELODIC --version"
+		$MELODIC --version
+		log_Debug_Msg "$MELODIC --help"
+		$MELODIC --help
+		log_Debug_Msg "md5sum $MELODIC"
+		md5sum $MELODIC
+	fi
+	log_Debug_Msg "End of melodic version log, help, and checksum"
+
+	# Use the Wishart-based dim estimate instead of the melodic-based estimate (--dim=${Dim})
+	# Also turn off melodic VN and feed in the input that has already been VN'ed instead (--vn flag, which turns OFF melodic's vn)
+	# Added "-m ${concatfmri}_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018
+	melodic_cmd=("${MELODIC}" -i "${concatfmrihp}_vnts" -o "${concatfmrihp}.ica/filtered_func_data.ica" --nobet --report --Oall --tr="${tr}" --vn --dim="${Dim}" -m "${concatfmri}_brain_mask")
+	if [ ! -z "${log_debugOn}" ] ; then
+		melodic_cmd+=(--verbose --debug)
+	fi
 
-#keep the existing handling code for melodic errors?
-debug_disable_trap
+	#keep the existing handling code for melodic errors?
+	debug_disable_trap
 
-log_Msg "melodic_cmd: ${melodic_cmd[*]}"
-"${melodic_cmd[@]}"
-return_code=$?
-log_Msg "melodic has been run: return_code = ${return_code}"
-log_Debug_Msg "melodic has been run: Contents of ${concatfmrihp}.ica follow"
-if [ ! -z "${log_debugOn}" ] ; then
-	ls -lRa ${concatfmrihp}.ica
-fi
+	log_Msg "melodic_cmd: ${melodic_cmd[*]}"
+	"${melodic_cmd[@]}"
+	return_code=$?
+	log_Msg "melodic has been run: return_code = ${return_code}"
+	log_Debug_Msg "melodic has been run: Contents of ${concatfmrihp}.ica follow"
+	if [ ! -z "${log_debugOn}" ] ; then
+		ls -lRa ${concatfmrihp}.ica
+	fi
 
-if [ "${return_code}" -ne "0" ] ; then
-	log_Err_Abort "melodic has returned a non-zero code"
-fi
+	if [ "${return_code}" -ne "0" ] ; then
+		log_Err_Abort "melodic has returned a non-zero code"
+	fi
 
-debug_enable_trap
+	debug_enable_trap
 
-# At this point (following melodic), ${concatfmrihp}_vnts is no longer necessary
-$FSLDIR/bin/imrm ${concatfmrihp}_vnts
+	# Delete some time series, resulting from the '--Oall' option in melodic, that aren't needed
+	# (these may only get created in the context of MIGP)
+	$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/alldat
+	$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/concat_data
 
-# Delete some time series, resulting from the '--Oall' option in melodic, that aren't needed
-# (these may only get created in the context of MIGP)
-$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/alldat
-$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/concat_data
+elif [[ $ICAmethod == ICASSO ]]; then
+	## ---------------------------------------------------------------------------
+	## Run icasso on concatenated file
+	## ---------------------------------------------------------------------------
+	
+	log_Msg "Running icasso ..."
+	matlab_code="${ML_PATHS} run_icasso('${Dim}','${concatfmri}','${concatfmrihp}','${ConcatFolder}','${icasso_vis}','${icasso_nICA}');exit"
+	log_Msg "Run interpreted MATLAB/Octave with code..."
+	log_Msg ${matlab_code}
+	
+	matlab -nodisplay -r "${matlab_code}";echo
 
+fi # if ICAmethod == MELODIC
+
+# At this point (following melodic or icasso), ${concatfmrihp}_vnts is no longer necessary
+$FSLDIR/bin/imrm ${concatfmrihp}_vnts
+	
 ## ---------------------------------------------------------------------------
 ## Housekeeping related to files expected for FIX
 ## ---------------------------------------------------------------------------
-
 cd "${concatfmrihp}.ica"
 
 # This is the concated volume time series from the 1st pass VN, with requested
@@ -829,7 +864,7 @@ fi
 
 # Other necessary files (N.B. Movement regressors already prep'ed above)
 $FSLDIR/bin/imln filtered_func_data.ica/mask mask
-if [ `$FSLDIR/bin/imtest ../${concatfmri}_SBRef` = 1 ] ; then
+if [ $($FSLDIR/bin/imtest ../${concatfmri}_SBRef) = 1 ] ; then
 	$FSLDIR/bin/imln ../${concatfmri}_SBRef mean_func
 else
 	$FSLDIR/bin/imln filtered_func_data.ica/mean mean_func
@@ -838,24 +873,21 @@ fi
 mkdir -p reg
 cd reg
 
-i_am_at=`pwd`
+i_am_at=$(pwd)
 log_Debug_Msg "current folder ${i_am_at}"
 
 $FSLDIR/bin/imln ../../../../T1w_restore_brain highres
 $FSLDIR/bin/imln ../../../../wmparc wmparc
 $FSLDIR/bin/imln ../mean_func example_func
 $FSLDIR/bin/makerot --theta=0 > highres2example_func.mat
-if [ `$FSLDIR/bin/imtest ../../../../T2w` = 1 ] ; then
+if [ $($FSLDIR/bin/imtest ../../../../T2w) = 1 ] ; then
 	$FSLDIR/bin/fslmaths ../../../../T1w -div ../../../../T2w veins -odt float
 	$FSLDIR/bin/flirt -in ../../../../brainmask_fs -ref veins -out veinbrainmask -applyxfm
 	$FSLDIR/bin/fslmaths veinbrainmask -bin veinbrainmask
-	$FSLDIR/bin/fslmaths veins -div `$FSLDIR/bin/fslstats veins -k veinbrainmask -P 50` -mul 2.18 -thr 10 -min 50 -div 50 veins
+	$FSLDIR/bin/fslmaths veins -div $($FSLDIR/bin/fslstats veins -k veinbrainmask -P 50) -mul 2.18 -thr 10 -min 50 -div 50 veins
 	$FSLDIR/bin/flirt -in veins -ref example_func -applyxfm -init highres2example_func.mat -out veins_exf
 	$FSLDIR/bin/fslmaths veins_exf -mas example_func veins_exf
 fi
-
-# Return to ${ConcatFolder}
-# Do not use 'cd ${ConcatFolder}', because ${ConcatFolder} may not be an absolute path
 cd ../..
 
 ## ---------------------------------------------------------------------------
@@ -863,7 +895,6 @@ cd ../..
 ## ---------------------------------------------------------------------------
 
 log_Msg "Running FIX"
-
 # Changes to handle user specified training data file
 if [[ "${TrainingData}" != "" ]]; then
 	# User has specified a training data file
@@ -905,6 +936,7 @@ else
 	# need to add a FIX version dependency to the script
 	fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}")
 fi
+
 log_Msg "fix_cmd: ${fix_cmd[*]}"
 ## MPH: The 'fix' script itself will continue to log to its own custom files
 ## Alert user to where those are
@@ -935,7 +967,7 @@ fi
 
 # The variance normalization ("_vn") outputs of fix (fix_3_clean) require use of fix1.067 or later
 # So check whether those files exist before moving/renaming them
-if [ `$FSLDIR/bin/imtest ${concatfmrihp}.ica/filtered_func_data_clean_vn` = 1 ]; then
+if [ $($FSLDIR/bin/imtest ${concatfmrihp}.ica/filtered_func_data_clean_vn) = 1 ]; then
 	$FSLDIR/bin/immv ${concatfmrihp}.ica/filtered_func_data_clean_vn ${concatfmrihp}_clean_vn
 fi
 if [ -f ${concatfmrihp}.ica/Atlas_clean_vn.dscalar.nii ]; then
@@ -984,8 +1016,8 @@ log_Msg "Splitting cifti and nifti back into individual runs"
 Start="1"
 for fmri in "${fmris[@]}" ; do
 	fmriNoExt=$($FSLDIR/bin/remove_ext $fmri)  # $fmriNoExt still includes leading directory components
-	NumTPS=`${Caret7_Command} -file-information ${fmriNoExt}_Atlas.dtseries.nii -no-map-info -only-number-of-maps`
-	Stop=`echo "${NumTPS} + ${Start} -1" | bc -l`
+	NumTPS=$(${Caret7_Command} -file-information ${fmriNoExt}_Atlas.dtseries.nii -no-map-info -only-number-of-maps)
+	Stop=$(echo "${NumTPS} + ${Start} -1" | bc -l)
 	log_Msg "${fmriNoExt}: Start=${Start} Stop=${Stop}"
 	
 	cifti_out=${fmriNoExt}_Atlas_hp${hp}_clean.dtseries.nii
@@ -1017,7 +1049,7 @@ for fmri in "${fmris[@]}" ; do
 		echo "$(tail -n ${nVols} ${fmriDir}/Movement_Regressors_hp${hp}_clean.txt)" > ${fmriDir}/Movement_Regressors_hp${hp}_clean.txt
 	fi
 
-	Start=`echo "${Start} + ${NumTPS}" | bc -l`
+	Start=$(echo "${Start} + ${NumTPS}" | bc -l)
 done
 
 cd "${DIR}"
diff --git a/ICAFIX/scripts/run_icasso.m b/ICAFIX/scripts/run_icasso.m
new file mode 100644
index 000000000..e9e81d171
--- /dev/null
+++ b/ICAFIX/scripts/run_icasso.m
@@ -0,0 +1,257 @@
+function sR = run_icasso(Dim,concatfmri,concatfmrihp,ConcatFolder,vis,nICA,maxIter)
+% run_icasso(Dim,concatfmri,concatfmrihp,ConcatFolder,vis,nICA,maxIter)
+% This function performs icasso decomposition for hcp_fix_multi_run.sh
+% It creates outputs in the style of MELODIC, so that FIX can be run
+% afterwards. Most but not all of MELODIC's outputs are created.
+% 
+% All of this function's inputs are strings
+%
+% Required Inputs (variable names verbatim from hcp_fix_multi_run.sh):
+%   Dim          : Data diminsionaliy estimate, typically Wishart-based
+%   concatfmri   : File name of original 4d time series without extension
+%   concatfmrihp : File name of high-passed 4d time series without extension
+%   ConcatFolder : Directory which contains concatfmri
+%
+% Optional Inputs:
+%   vis     : Whether to create and save icasso figures, see icasso.m 'basic' (default) or 'off'
+%   nICA    : Number of ICA repetitions per icasso repetition, @ delimited string, (default = '100') 
+%             The number of icasso repetitions is equal to the the number delimiters + 1
+%   maxIter : Maximum number of iterations per ica fit (default = '1000')
+%
+% example inputs:
+% Dim = '41';
+% concatfmri = 'tfMRI_EMOTION_RL_LR';
+% concatfmrihp = 'tfMRI_EMOTION_RL_LR_hp0';
+% ConcatFolder = '/mnt/myelin/burke/HCPpipelines/dev_study/100307/MNINonLinear/Results/tfMRI_EMOTION_RL_LR_ICASSO';
+% vis = 'basic';
+% nICA = '2@2';
+% maxIter = '1000';
+% 
+% Dim = str2double(Dim);
+% nICA = cellfun(@str2double,regexp(nICA,'@','split'));
+% maxIter = str2double(maxIter);
+
+% Created 2024-06-05
+% Burke Rosen
+
+% ToDo:
+% () Currently the vnts file and brainMaskFile paths and names are inferred from
+%    ConcatFolder, concatfmrihp, and concatfmri. Maybe they should be their own 
+%    arguments for flexibility.
+
+%% parse parameters
+if nargin < 4 || any(cellfun(@isempty,{Dim,concatfmri,concatfmrihp,ConcatFolder}))
+  error('Dim, concatfmri, concatfmrihp, and ConcatFolder are required!')
+end
+Dim = str2double(Dim);
+if nargin < 5 || isempty(vis) 
+  vis = 'basic';
+end
+if ~ismember(vis,{'basic','off'})
+  warning('vis must be ''basic'' or ''off'', reverting to ''basic''')
+  vis = 'basic';
+end
+if nargin < 6 || isempty(nICA) 
+  nICA = 100;
+else
+  nICA = cellfun(@str2double,regexp(nICA,'@','split'));
+end
+if nargin < 7 || isempty(maxIter) 
+  maxIter = 1000;
+else
+  maxIter = str2double(maxIter);
+end
+
+%% check IO dependencies
+% use imaging processing toolbox utilities, or if those are not available use FSL utilities 
+function out = out2(fun);[~,out] = fun();end
+function out = out3(fun);[~,~,out] = fun();end
+if isempty(which('niftiread'))
+  if isempty(which('read_avw'))
+    error('neither niftiread nor read_avw on matlab path!')
+  end
+  infoNIFTI = @(fName) struct('ImageSize',out2(@() read_avw(fName))','PixelDimensions',out3(@() read_avw(fName))');
+  readNIFTI = @(fName) read_avw(fName);
+  writeNIFTI = @(img,fName,hdr) save_avw(img,fName,'f',hdr.PixelDimensions);
+else
+  readNIFTI = @(fName) niftiread(fName);
+  infoNIFTI = @(fName) niftiinfo(fName);
+  writeNIFTI = @(img,fName,hdr) niftiwrite(img,fName,hdr,'Compressed',true);
+end
+
+%% parse paths
+% inputs
+vntsFile = sprintf('%s/%s_vnts.nii.gz',ConcatFolder,concatfmrihp);
+brainMaskFile = sprintf('%s/%s_brain_mask.nii.gz',ConcatFolder,concatfmri);
+if ~exist(vntsFile,'file');error('%s doesn''t exist!',vntsFile);end
+if ~exist(brainMaskFile,'file');error('%s doesn''t exist!',brainMaskFile);end
+
+% outputs
+outDir = sprintf('%s/%s.ica/filtered_func_data.ica',ConcatFolder,concatfmrihp);
+if ~mkdir(outDir);error('Unable to make output folder!');end
+
+%% load data, reshape to 2d, and apply brain mask
+vnts = double(readNIFTI(vntsFile));% fastica needs double
+brainMask = logical(readNIFTI(brainMaskFile));
+volDim = size(vnts);
+[vnts,mtxDim] = maskAndSpatiallyFlatten(vnts,brainMask);
+
+%% Wishart filter
+% nWF = 1;
+% icaDimOut = icaDim(vnts,-1,1,-1,nWF);
+% vnts = double(icaDimOut.data);
+
+%% remove mean image and check for voxels with constant values
+% vntsT = vnts';
+% vntsT = vntsT - mean(vntsT);
+% if any(all(vntsT == vntsT(1,:)))
+%   error('some voxels have constant value!')
+% end
+
+%% run icasso
+% note: icasso fixes the randomization seed with rng('default')
+sR = cell(numel(nICA),1);
+[iq,A,~,~,sR{1}] = ...
+  icasso('both',vntsT,nICA(1),'approach','symm','g','pow3',...
+  'lastEig',Dim,'numOfIC',Dim,'maxNumIterations',maxIter,'vis',vis); 
+if strcmp(vis,'basic'); printFigs(outDir,1);end
+for iC = 2:numel(nICA)
+  [iq,A,~,~,sR{iC}] = ...
+    icasso('bootstrap',vntsT,nICA(iC),'approach','symm','g','pow3','initGuess',A,...
+    'lastEig',Dim,'numOfIC',Dim,'maxNumIterations',maxIter,'vis',vis); 
+  if strcmp(vis,'basic'); printFigs(outDir,iC);end
+end
+[pcaE,pcaD] = fastica(vntsT,'only','pca');
+[S_final,A_final,~] = ...
+  fastica(vntsT,'initGuess',A,'approach','symm','g','pow3',...
+  'lastEig',Dim,'numOfIC',Dim,'pcaE',pcaE,'pcaD',pcaD,...
+  'displayMode','off','maxNumIterations',maxIter);
+pcaD = diag(pcaD);
+totVariance = sum(pcaD(end-Dim+1:end))./sum(pcaD);
+  % proportion of total variance explained by first <Dim> components
+clear vntsT
+S_final = S_final';
+% combineMultiLevelIcassoFigs
+% keyboard
+%% set component sign and sort by variance explained
+% set sign
+maxSfinal = max(S_final);
+absminSfinal = abs(min(S_final));
+pos = maxSfinal > absminSfinal;
+neg = maxSfinal < absminSfinal;
+signAll = sign(pos + neg * -1);
+S_final = S_final .* signAll;
+A_final = A_final .* signAll;
+
+% sort by variance explained
+tVar = var(A_final);
+tICAPercentVariances = tVar / sum(tVar) * 100;
+[tICAPercentVariances,Is] = sort(tICAPercentVariances,'descend');
+S_final = S_final(:,Is);
+A_final = A_final(:,Is);
+iq = iq(Is);
+
+W_final = pinv(A_final);
+
+%% calculate single-regression z-stats
+NODEtsnorm = normalise(A_final);% z-score
+pN = pinv(NODEtsnorm); dpN = diag(pN * pN')';
+dof = size(NODEtsnorm,1) - size(NODEtsnorm,2) - 1;
+residuals = demean(vnts, 2) - S_final * NODEtsnorm';
+t = double(S_final ./ sqrt(sum(residuals .^ 2, 2) * dpN / dof));
+Z = zeros(size(t));
+Z(t > 0) = min(-norminv(tcdf(-t(t > 0), dof)),  38.5);
+Z(t < 0) = max( norminv(tcdf( t(t < 0), dof)), -38.5);
+Z(isnan(Z)) = 0;
+
+%% save icasso outputs to match melodic
+% reshape betas (S_final) and z-scores back to 4-d
+mtxDim = [mtxDim(1) Dim];
+volDim = [volDim(1:3) Dim];
+S_final4d = unmaskAndSpatiallyInflate(S_final,volDim,brainMask,mtxDim);
+Z4d = unmaskAndSpatiallyInflate(Z,volDim,brainMask,mtxDim);
+
+% save betas and z-score as 4-d nifti volumes
+hdr = infoNIFTI(vntsFile);
+hdr.ImageSize(end) = Dim;
+writeNIFTI(single(S_final4d),sprintf('%s/melodic_oIC',outDir),hdr)
+writeNIFTI(single(Z4d),sprintf('%s/melodic_IC',outDir),hdr)
+
+% save mixing matrices as tab-delimited text
+dlmwrite(sprintf('%s/melodic_unmix', outDir), W_final, '\t');
+dlmwrite(sprintf('%s/melodic_mix',   outDir), A_final, '\t');
+copyfile(sprintf('%s/melodic_mix',   outDir), sprintf('%s/melodic_Tmodes', outDir));% mix and Tmodes are the same
+
+% save variance explained as tab-delimited text
+ICstats = [tICAPercentVariances tICAPercentVariances * totVariance]; 
+dlmwrite([outDir '/melodic_ICstats'], ICstats, 'delimiter', '\t');
+
+% copy brainmask into melodic dir
+copyfile(brainMaskFile,[outDir '/mask.nii.gz'])
+
+% save pcaD and pcaE? Looks like FIX doesn't need them
+
+%% calculate and save FTmix spectra
+ts = struct();
+ts.Nsubjects = 1;
+ts.ts = A_final;
+[ts.NtimepointsPerSubject,ts.Nnodes] = size(ts.ts);
+ts_spectra = nets_spectra_sp(ts);
+dlmwrite([outDir '/melodic_FTmix'], ts_spectra, 'delimiter', '\t');
+fid = fopen([outDir '/components.txt'],'w');
+fprintf(fid,'%i: Signal\n',1:size(S_final,2));fclose(fid);
+
+% FTmix.sdseries.nii aren't need by FIX, so don't produce them
+% [~,~] = system(['wb_command -cifti-create-scalar-series ' ...
+%   sprintf('%s/melodic_FTmix %s/melodic_FTmix.sdseries.nii -transpose -name-file %s/components.txt -series HERTZ 0 %f',...
+%   outDir,outDir,outDir,tr)]);
+
+%% Mixture modeling
+fprintf('performing melodic mixture modeling ...\n')
+mixtureModel([outDir '/melodic_IC.nii.gz'],[outDir '/melodic_IC.nii.gz']);% overwrites in-place without saving full report
+
+%% Helper subfunctions
+function printFigs(outD,lvl)
+  figH = findall(0,'type','figure');
+  figH = figH(~contains({figH.Name},'centrotypes'));% the centrotypes figure isn't useful
+  for iF = 1:numel(figH)
+    figFile = sprintf('%s/%s_%i.fig',...
+      outD,strrep(strrep(strrep(figH(iF).Name,' ','_'),':',''),'Icasso','icasso'),lvl);
+    try
+      savefig(figH(iF),figFile,'compact');
+    catch 
+      warning('Could not save icasso figures.')
+    end
+  end
+  close all;
+end
+
+function [mtx,mtxDim] = maskAndSpatiallyFlatten(img,msk)
+  % reshapes 4d img into 2d matrix where the first 3 dims are put into the 1st dim
+  % second arg msk is 3d logical volume
+  % second output is size before mask is applied
+
+  imgDim = size(img);
+  mtx = reshape(img,prod(imgDim(1:3)),imgDim(4));
+  mtxDim = size(mtx);%
+  if nargin == 2
+    mtx = mtx(logical(msk(:)),:);
+  end
+end
+
+function img = unmaskAndSpatiallyInflate(mtx,imgDim,msk,mtxDim)
+  % reshapes 2d mtx into 4d matrix where the first 1 dim are put into the
+  % 1st 3 dims using the dimensions supplied in imgDim
+  % if a msk is supplied then the size of the 2d mtx before the mask was
+  % applied is needed
+
+  if nargin > 2
+    img = zeros(mtxDim);
+    img(msk,:) = mtx;
+  else
+    img = mtx;
+  end
+  img = reshape(img,imgDim);
+end
+
+end %EOF
\ No newline at end of file
diff --git a/fMRIVolume/scripts/multiEchoCombine.m b/fMRIVolume/scripts/multiEchoCombine.m
index a19f558c8..b4c9e55f4 100755
--- a/fMRIVolume/scripts/multiEchoCombine.m
+++ b/fMRIVolume/scripts/multiEchoCombine.m
@@ -53,7 +53,8 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 % Created 2024-03-19
 % Burke Rosen 
 
-%% handle inputs, set defaults, and load data
+
+%% handle inputs, set defaults
 if ~exist(tscPath,'file')
   error('%s does not exist!',tscPath)
 end
@@ -94,11 +95,28 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
   nonlinAlgo = 'Levenberg-Marquardt'; 
 end
 
-% load data
-I = niftiread(tscPath);
-hdr = niftiinfo(tscPath);
+%% set nifti IO functions
+% use imaging processing toolbox utilities, or if those are not available, use FSL utilities 
+function out = out2(fun);[~,out] = fun();end
+function out = out3(fun);[~,~,out] = fun();end
+if isempty(which('niftiread'))
+  if isempty(which('read_avw'))
+    error('neither niftiread nor read_avw on matlab path!')
+  end
+  infoNIFTI = @(fName) struct('ImageSize',out2(@() read_avw(fName))','PixelDimensions',out3(@() read_avw(fName))');
+  readNIFTI = @(fName) read_avw(fName);
+  writeNIFTI = @(img,fName,hdr) save_avw(img,fName,'f',hdr.PixelDimensions);
+else
+  readNIFTI = @(fName) niftiread(fName);
+  infoNIFTI = @(fName) niftiinfo(fName);
+  writeNIFTI = @(img,fName,hdr) niftiwrite(img,fName,hdr,'Compression',true);
+end
+
+%% load data
+I = readNIFTI(tscPath);
+hdr = infoNIFTI(tscPath);
 if ~isempty(sctPath)
-  S = niftiread(sctPath);
+  S = readNIFTI(sctPath);
 end
 
 % get dims
@@ -118,7 +136,7 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
     Y(:,:,:,iE) = mean(I(:,:,:,(1:framePerEcho)+framePerEcho*(iE-1)),4);
   end
 else
-  Y = niftiread(refPath);
+  Y = readNIFTI(refPath);
 end
 
 %% Create time matrix
@@ -265,8 +283,8 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 hdr3d.PixelDimensions = hdr3d.PixelDimensions(1:3);
 T2starPath = strrep(tscPath,'.nii.gz','_T2star');
 S0Path = strrep(tscPath,'.nii.gz','_S0');
-niftiwrite(T2star,T2starPath,hdr3d,'Compressed',true)
-niftiwrite(S0,S0Path,hdr3d,'Compressed',true)
+writeNIFTI(T2star,T2starPath,hdr3d)
+writeNIFTI(S0,S0Path,hdr3d)
 
 %% Extrapolate voxels with negative T2* from T2* and S0
 % and voxels with all echoes below noisefloor, if fitNoiseFloor
@@ -274,7 +292,7 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 badVoxel = single(T2star < 0);
 if ~isempty(find(badVoxel,1))% any(x,'all') syntax isn't present in 2017b
   bvPath = tempname;
-  niftiwrite(badVoxel,bvPath,hdr3d,'Compressed',true)
+  writeNIFTI(badVoxel,bvPath,hdr3d)
   cmd = sprintf(...
     'wb_command -volume-dilate %s.nii.gz 5 WEIGHTED %s.nii.gz -bad-voxel-roi %s.nii.gz',...-grad-extrapolate
     T2starPath,T2starPath,bvPath);
@@ -283,7 +301,7 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
     'wb_command -volume-dilate %s.nii.gz 5 WEIGHTED %s.nii.gz -bad-voxel-roi %s.nii.gz',...-grad-extrapolate
     S0Path,S0Path,bvPath);
   [~,~] = unix(cmd);
-  T2star = niftiread(T2starPath);
+  T2star = readNIFTI(T2starPath);
   delete([bvPath '.nii.gz'])
 end
 clear S0;
@@ -307,7 +325,7 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 hdrW = hdr;
 hdrW.ImageSize(4) = nE;
 WeightsPath = strrep(tscPath,'.nii.gz','_EchoWeights');
-niftiwrite(W,WeightsPath,hdrW,'Compressed',true)
+writeNIFTI(W,WeightsPath,hdrW)
 
 %% apply weights and export results
 if ~isempty(sctPath)
@@ -317,7 +335,7 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 
   % save combined Scout / SBref
   hdr.ImageSize(4) = framePerEcho;
-  niftiwrite(S,strrep(sctPath,'.nii.gz','_CombEchoes'),hdr3d,'Compressed',true)
+  writeNIFTI(S,strrep(sctPath,'.nii.gz','_CombEchoes'),hdr3d)
 end
 
 % apply to timeseries
@@ -327,6 +345,6 @@ function multiEchoCombine(tscPath,TE,refPath,sctPath,T2starMethod,weightMethod,f
 
 % save combined timeseries
 hdr.ImageSize(4) = framePerEcho;
-niftiwrite(I,strrep(tscPath,'.nii.gz','_CombEchoes'),hdr,'Compressed',true)
+writeNIFTI(I,strrep(tscPath,'.nii.gz','_CombEchoes'),hdr)
 
 end % EOF
diff --git a/global/matlab/mixtureModel.m b/global/matlab/mixtureModel.m
new file mode 100644
index 000000000..40330d627
--- /dev/null
+++ b/global/matlab/mixtureModel.m
@@ -0,0 +1,93 @@
+function mixtureModel(inFile,outFile,wbcmd,melocmd)
+% mixtureModel(inFile,outFile,wbcmd,melocmd)
+% Performs Gaussian mixture modeling on precomputed ICs
+% Wrapper arround melodic.Accepts nifti or cifti inputs. 
+% Require Input:
+%   inFile  : file path to IC z-scores, file path including extension as string (accepts nifti or cifti files)
+%   outFile : file path to IC z-scores with Gaussian mixture modeling
+%              Can output nifti, nifit-gz, or cifti, depending on extension: .nii, .nii.gz, .dscalar.nii,
+%              Output format does not have to match input format, unless output is cifti.
+% Optional Inputs
+%   wbcmd   : worbench command, defaults to 'wb_command' 
+%   melocmd : melodic command, defaults to 'melodic'
+%
+% See also: 
+% https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=fsl;6e85d498.1607
+% https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/MELODIC#Using_melodic_for_just_doing_mixture-modelling
+% https://www.fmrib.ox.ac.uk/datasets/techrep/tr02cb1/tr02cb1.pdf
+
+% Created by Burke Rosen
+% 2024-07-08
+% Dependencies:workbench, FSL
+% Written with workbench 2.0 and FSL 6.0.7.1
+%
+% ToDo: 
+% Feed arbitrary melodic arguments. 
+
+%% handle inputs 
+if nargin < 2 || isempty(inFile) || isempty(outFile)
+  error('inFile and outFile arguments are required!'); 
+end 
+if nargin < 3 || isempty(wbcmd); wbcmd = 'wb_command'; end 
+if nargin < 4 || isempty(melocmd); melocmd = 'melodic'; end 
+[wbStat,~] = system(wbcmd);
+[meloStat,~] = system(['which ' melocmd]);
+if wbStat; error('workbench_command binary %s not on path',wbcmd);end
+if meloStat; error('melodic command binary %s not on path',melocmd);end
+inFile0 = inFile;
+outFile0 = outFile;
+tDir = tempname;
+FSLOUTPUTTYPE0 = getenv('FSLOUTPUTTYPE');
+stat = zeros(9,1);
+if endsWith(inFile0,'dscalar.nii')
+  % convert to input cifti to nifti
+  inFile = [tempname '.nii'];
+  [stat(1),~] = system(sprintf('%s -cifti-convert -to-nifti %s %s -smaller-dims',wbcmd,inFile0,inFile),'-echo');
+  [stat(2),out] = system(sprintf('%s -file-information %s | grep Dimensions',wbcmd,inFile),'-echo');% more elegant to use niftiinfo, but that adds dependency
+  dims = str2num(out(12:end));
+  if any(dims == 1) && find( dims == 1,1) < 4
+    warning('singleton dimension in converted nifti, melodic may not interpret correctly!');
+  end
+  inFile = strrep(inFile,'.nii','');
+elseif endsWith(inFile0,'.nii') 
+  inFile = strrep(inFile0,'.nii','');
+elseif endsWith(inFile0,'.nii.gz') 
+  inFile = strrep(inFile0,'.nii.gz','');
+else
+  error('inFile is not a nifti or dscalar cifti?')
+end
+if endsWith(outFile0,'dscalar.nii')
+  if ~endsWith(inFile0,'dscalar.nii')
+    error('cifti output only supported for cifti input.')
+  end
+  outFile = strrep(outFile0,'.dscalar.nii','');
+  FSLOUTPUTTYPE = 'NIFTI_GZ';
+elseif endsWith(outFile0,'.nii') 
+  outFile = strrep(outFile0,'.nii','');
+  FSLOUTPUTTYPE = 'NIFTI';
+elseif endsWith(outFile0,'.nii.gz') 
+  FSLOUTPUTTYPE = 'NIFTI_GZ';
+  outFile = strrep(outFile0,'.nii.gz','');
+else
+  error('outFile is not a nifti or dscalar cifti?')
+end
+
+%% run gaussian mixture modeling with melodic
+[stat(3),~] = system(sprintf('mkdir -p %s;echo "1" > %s/grot', tDir, tDir),'-echo');
+[stat(4),~] = system(sprintf(...
+  '%s -i %s --ICs=%s --mix=%s/grot -o %s --Oall --report -v --mmthresh=0',... 
+  melocmd,inFile, inFile, tDir, tDir),'-echo');
+[stat(5),~] = system(sprintf(...
+    'FSLOUTPUTTYPE=%s;fslmerge -t %s $(ls %s/stats/thresh_zstat* | sort -V);FSLOUTPUTTYPE=%s;',...
+    FSLOUTPUTTYPE, outFile, tDir ,FSLOUTPUTTYPE0),'-echo');
+[stat(6),~] = system(['rm -r ' tDir],'-echo');% clean up temporary files
+
+%% convert output to cifti, if needed
+if endsWith(outFile0,'dscalar.nii')
+  [stat(7),~] = system(sprintf('%s -cifti-convert -from-nifti %s.nii.gz %s %s',wbcmd,outFile,inFile0,outFile0),'-echo');
+  [stat(8),~] = system(sprintf('imrm %s %s',inFile,outFile),'-echo');% clean up intermediate niftis
+elseif endsWith(inFile0,'dscalar.nii')
+  [stat(9),~] = system(sprintf('imrm %s',inFile),'-echo');
+end
+if any(stat);error('system commands failed, see console output!');end
+end % EOF
\ No newline at end of file