function DEM = DEM_cells

% This demo illustrates self organisation in an ensemble of (sixteen) cells
% using the same principles described in DEM_morphogenesis, using a simpler
% generative model. Overall, the dynamics of these simulations show how one
% can prescribe a point attractor for each constituent of an ensemble that
% endows the ensemble with a point attractor to which the ensemble
% converges. In this example, we consider the special case where the point
% attractor is itself a Markov blanket. In other words, cells come to
% acquire dependencies, in terms of intracellular signalling, that conform
% to a simple Markov blanket with intrinsic or internal cells, surrounded
% by active cells that are, in turn, surrounded by sensory cells. This
% organisation rests upon intracellular signals and active inference using
% generalised (second-order) variational filtering. In brief, the hidden
% causes driving action (migration and signalling) are expectations about
% cell type. These expectations are optimised using sensory signals;
% namely, the signals generated by other cells. By equipping each cell with
% prior beliefs about what it would sense if it was a particular cell type
% (i.e., internal, active or sensory), the act (i.e., move and signal) so
% that they behave and infer their role in an ensemble of cells that itself
% has a Markov blanket. In a DEM_cell_cell.m, we use this first-order
% scheme to simulate hierarchical emergence of Markov blankets; i.e.,
% ensembles of cells that can be one of three types at the local level;
% independently of their role at the global level.
%__________________________________________________________________________
% Copyright (C) 2017 Wellcome Trust Centre for Neuroimaging
 
% Karl Friston
% $Id: DEM_cells.m 7150 2017-08-08 19:57:33Z karl $

% preliminaries
%--------------------------------------------------------------------------
clear global
rng('default')
N         = 32;                              % duration of process

% generative process and model
%==========================================================================
M(1).E.d  = 1;                               % approximation order
M(1).E.n  = 2;                               % embedding order
M(1).E.s  = 1;                               % smoothness
M(1).E.dt = 2;                               % time step

% Target morphology
%==========================================================================
n     = [1 6 9];                             % internal, active and sensory
r     = [0 1 3];                             % radial distance from centre
pos   = [];                                  % position
typ   = [];                                  % cell type
for j = 1:length(n)
    for i = 1:n(j)
        pos(:,end + 1) = r(j)*[sin(2*pi*i/n(j));cos(2*pi*i/n(j))];
        typ(j,end + 1) = 1;
    end
end

% Target sensation - given the canonical location and types
%--------------------------------------------------------------------------
sen   = foxhound(pos,typ);
for i = 1:3
    P.sen(:,i) = mean(sen(:,find(typ(i,:))),2);
end

% check the predicted sensations are sufficiently orthogonal for inference
%--------------------------------------------------------------------------
disp(P.sen)

% initialise expectations and action
%--------------------------------------------------------------------------
v     = typ/2 + randn(size(typ))/8;          % states (identity)                          % predicted sensations
a.pos = pos/2 + randn(size(pos))/8;          % chemotaxis
a.sec = spm_softmax(v);                      % secretion
g     = Mg([],v,P);
   
n     = sum(n);                              % number of cells
ns(1) = size(g.sec,1);                       % number of secretions
ns(2) = size(g.sen,1);                       % number of sensations
na(1) = size(a.pos,1);                       % number of positions
na(2) = size(a.sec,1);                       % number of secretions

% restriction matrix (ensuring action is strictly local)
%--------------------------------------------------------------------------
In    = eye(n,n);
R     = spm_cat({kron(In,ones(ns(1),na(1))) kron(In,ones(ns(1),na(2)));
                 kron(In,ones(ns(2),na(1))) kron(In,ones(ns(2),na(2)))});

% level 1 of generative process
%--------------------------------------------------------------------------
G(1).g  = @(x,v,a,P) Gg(x,v,a,P);
G(1).v  = Gg([],[],a,a);
G(1).V  = exp(16);                           % precision (noise)
G(1).U  = exp(-2);                           % precision (action)
G(1).R  = R;                                 % rate matrix
G(1).pE = a;                                 % form (action)
G(1).aP = exp(-8);                           % precision (action)


% level 2; causes (action)
%--------------------------------------------------------------------------
G(2).a  = spm_vec(a);                        % endogenous cause (action)
G(2).v  = 0;                                 % exogenous  cause
G(2).V  = exp(16);


% generative model
%==========================================================================

% level 1 of the generative model:
%--------------------------------------------------------------------------
M(1).g  = @(x,v,P) Mg([],v,P);
M(1).v  = g;
M(1).V  = exp(4);                            % precision of sensations
M(1).pE = P;

% level 2:
%--------------------------------------------------------------------------
M(2).v  = v;
M(2).V  = exp(-4);                           % prior precision of identity

% hidden cause and prior identity expectations (and time) - none
%--------------------------------------------------------------------------
U     = zeros(n*3,N);
C     = zeros(1,N);

% assemble model structure
%--------------------------------------------------------------------------
DEM.M = M;
DEM.G = G;
DEM.C = C;
DEM.U = U;

% solve
%==========================================================================
DEM   = spm_ADEM(DEM);
spm_DEM_qU(DEM.qU,DEM.pU);


% Graphics
%==========================================================================

% Evolution
% --------------------------------------------------------------------------
subplot(2,2,2)
for t = 1:N
    v = spm_unvec(DEM.qU.a{2}(:,t),a);
    c = spm_unvec(DEM.qU.v{2}(:,t),DEM.M(2).v);
    for i = 1:n
        x = v.pos(1,i);
        y = v.pos(2,i);
        R = spm_softmax(c(:,i),2);
        plot(x,y,'o','markersize',4,'markeredgecolor',[0 0 0],'markerfacecolor',R); hold on
        if t == N
            plot(x,y,'o','markersize',10,'markeredgecolor',[0 0 0],'markerfacecolor',R); hold on
        end
    end
    title ('Self-organisation','fontsize',15);
    xlabel('position');
    ylabel('position');
    axis([-1 1 -1 1]*r(end))
    axis equal off, drawnow
end

% subroutines
%==========================================================================

% Generating sensations: self signalling and extracellular signalling
%--------------------------------------------------------------------------
function g = Gg(x,v,a,P)

a      = spm_unvec(a,P);                   % action
g.sec  = a.sec;                            % secretion
g.sen  = foxhound(a.pos,a.sec);            % sensation

return

% Generating predictions based on expectations about subtype
%--------------------------------------------------------------------------
function g = Mg(x,v,P)

p      = spm_softmax(v);                   % expected identity
g.sec  = p;                                % secretion
g.sen  = P.sen*p;                          % sensation

return

% Sensed signals that decay as a Gaussian function of distance
%--------------------------------------------------------------------------
function sen = foxhound(x,y)

% sen = sensation
% x   = position
% y   = cell or signal type (internal, active sensory)

[m,n] = size(y);                          % number of signals and cells
sen   = zeros(m,n);                       % extracellular signals
k     = 2;                                % spatial decay constant
for i = 1:n
    for j = 1:n
        
        % distance
        %------------------------------------------------------------------
        d = x(:,i) - x(:,j);
        d = d'*d;
        
        % signal sensed
        %------------------------------------------------------------------
        if i ~= j
           sen(:,i) = sen(:,i) + exp(-k*d)*y(:,j);
        end

    end
end
return