function plotSummary ( D )
%plots summary statistics for the repetitions
%D = zeros ( NReps, MaxTime+1, Npops );

[ nReps maxTime nPops ] = size ( D );

%extract the maximum population size (all repetitions, all time points, all
%populations) for axis scaling
TimePts = 0 : maxTime-1;

%calculate the mean and variance between populations for each repetition
MeanP = zeros ( nReps, maxTime ); %mean population size
VarP = zeros ( nReps, maxTime ); %variance in population sizes
PrevP = zeros ( nReps, maxTime ); %proportion of populations that are non-zero (not extinct)
MPSize = zeros ( nReps, maxTime ); %metapopulation size
for irep = 1 : nReps
    Dr = squeeze ( D ( irep, :, : ) ); %Dr is a 2-D array of times and populations
    MeanP ( irep, : ) = mean ( Dr, 2 ); %mean over time
    VarP ( irep, : ) =  var ( Dr, 0, 2 ); %variance over time
    MPSize ( irep, : ) = sum ( Dr, 2 ); %sum of population sizes
    for itime = 1 : maxTime
        index = find ( Dr ( itime, : ) ); %index to the non-zero elements
        PrevP ( irep, itime ) = length ( index ) / nPops;
    end
end

%plots summary statistics for all simulations
figure
subplot ( 2, 2, 1 )
plot ( TimePts, MeanP' ) %mean of each simulation
hold on
plot ( TimePts, mean ( MeanP ), 'kx' ) %mean of the means
xlabel ( 'Time' )
ylabel ( 'Mean Pop. Size' )

subplot ( 2, 2, 2 )
plot ( TimePts, MPSize' ); %metapopulation size
hold on
plot ( TimePts, mean ( MPSize ), 'kx' )
xlabel ( 'Time' );
ylabel ( 'Metapop. Size' );

subplot ( 2, 2, 3 )
plot ( TimePts, VarP' ) %variance of each simulation
hold on
plot ( TimePts, mean ( VarP ), 'kx' ) %mean of the variances
xlabel ( 'Time' )
ylabel ( 'Variance' )
subplot ( 2, 2, 3 )

subplot ( 2, 2, 4 )
plot ( TimePts, PrevP' ) %proportion of populations surviving
hold on
plot ( TimePts, mean ( PrevP ), 'kx' ) %mean of surviving proportion
axis ( [ 0 maxTime-1 0 1 ] );
xlabel ( 'Time' )
ylabel ( 'Prop Pops Surviving' );