Modeling and MEG evidence of early consonance processing in auditory cortex (Tabas et al 2019)

 Download zip file 
Help downloading and running models
Accession:256624
Pitch is a fundamental attribute of auditory perception. The interaction of concurrent pitches gives rise to a sensation that can be characterized by its degree of consonance or dissonance. In this work, we propose that human auditory cortex (AC) processes pitch and consonance through a common neural network mechanism operating at early cortical levels. First, we developed a new model of neural ensembles incorporating realistic neuronal and synaptic parameters to assess pitch processing mechanisms at early stages of AC. Next, we designed a magnetoencephalography (MEG) experiment to measure the neuromagnetic activity evoked by dyads with varying degrees of consonance or dissonance. MEG results show that dissonant dyads evoke a pitch onset response (POR) with a latency up to 36 ms longer than consonant dyads. Additionally, we used the model to predict the processing time of concurrent pitches; here, consonant pitch combinations were decoded faster than dissonant combinations, in line with the experimental observations. Specifically, we found a striking match between the predicted and the observed latency of the POR as elicited by the dyads. These novel results suggest that consonance processing starts early in human auditory cortex and may share the network mechanisms that are responsible for (single) pitch processing.
Reference:
1 . Tabas A, Andermann M, Schuberth V, Riedel H, Balaguer-Ballester E, Rupp A (2019) Modeling and MEG evidence of early consonance processing in auditory cortex. PLoS Comput Biol 15:e1006820 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Auditory cortex;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB; Python;
Model Concept(s): Magnetoencephalography;
Implementer(s): Tabas, Alejandro [tabas at cbs.mpg.de];
function s = tdochCortex(r, pars)

    if nargin == 1, pars = loadParameters(); end
  
    s = initialiseStates(r, pars);
    s = corticalProcessing(r, s, pars);

end



function s = corticalProcessing(r, s, pars)

   if pars.verb, fprintf('Simulating cortical processing...');  end;

    nOfIts = length(r.timeSpace) - 1;
    markerThresh = nOfIts / 20;
    marker = 0;

    for ti = 1:(length(r.timeSpace) - 1)

        if pars.verb
            marker = marker + 1;
            if marker > markerThresh
                fprintf('.');
                marker = 0;
            end
        end

        dt = r.timeSpace(ti + 1) - r.timeSpace(ti);
        s = updateStateVars(s, r, dt, pars);

    end

    if pars.verb, fprintf('..done! <3 time = %.1fm\n', toc / 60); end;

end



function s = initialiseStates(r, pars)

    % Parameters
    outOfClockIterations = 2000;

    dt = r.timeSpace(3) - r.timeSpace(2);

    s0.t = 1;
    s0.n.SPn = zeros(outOfClockIterations + 1, pars.N);
    s0.n.SPa = zeros(outOfClockIterations + 1, pars.N);

    s0.p.Sg  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.Sn  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.Sa  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.Hi  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.He  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.xi  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.xe  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.Ai  = zeros(outOfClockIterations + 1, pars.N);
    s0.p.Ae  = zeros(outOfClockIterations + 1, pars.N);
    
    s0.q.Sg  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.Sn  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.Sa  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.Hi  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.He  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.xi  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.xe  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.Ai  = zeros(outOfClockIterations + 1, pars.N);
    s0.q.Ae  = zeros(outOfClockIterations + 1, pars.N);

    r0 = r;
    r0.A = zeros(outOfClockIterations + 1, pars.N);
    r0.E = zeros(outOfClockIterations + 1, 1);

    for ii = 1:outOfClockIterations
        s0 = updateStateVars(s0, r0, dt, pars);
    end

    s.n.SPn = zeros(length(r.timeSpace), pars.N);
    s.n.SPa = zeros(length(r.timeSpace), pars.N);

    s.p.Sg  = zeros(length(r.timeSpace), pars.N);
    s.p.Sn  = zeros(length(r.timeSpace), pars.N);
    s.p.Sa  = zeros(length(r.timeSpace), pars.N);
    s.p.Hi  = zeros(length(r.timeSpace), pars.N);
    s.p.He  = zeros(length(r.timeSpace), pars.N);
    s.p.xi  = zeros(length(r.timeSpace), pars.N);
    s.p.xe  = zeros(length(r.timeSpace), pars.N);
    s.p.Ai  = zeros(length(r.timeSpace), pars.N);
    s.p.Ae  = zeros(length(r.timeSpace), pars.N);

    s.q.Sg  = zeros(length(r.timeSpace), pars.N);
    s.q.Sn  = zeros(length(r.timeSpace), pars.N);
    s.q.Sa  = zeros(length(r.timeSpace), pars.N);
    s.q.Hi  = zeros(length(r.timeSpace), pars.N);
    s.q.He  = zeros(length(r.timeSpace), pars.N);
    s.q.xi  = zeros(length(r.timeSpace), pars.N);
    s.q.xe  = zeros(length(r.timeSpace), pars.N);
    s.q.Ai  = zeros(length(r.timeSpace), pars.N);
    s.q.Ae  = zeros(length(r.timeSpace), pars.N);

    s.t = 1;
    s.n.SPn(1, :) = s0.n.SPn(end, :); 
    s.n.SPa(1, :) = s0.n.SPa(end, :); 

    s.p.Sg(1, :)  = s0.p.Sg(end, :); 
    s.p.Sn(1, :)  = s0.p.Sn(end, :); 
    s.p.Sa(1, :)  = s0.p.Sa(end, :); 
    s.p.Hi(1, :)  = s0.p.Hi(end, :);
    s.p.He(1, :)  = s0.p.He(end, :);
    s.p.xi(1, :)  = s0.p.xi(end, :);
    s.p.xe(1, :)  = s0.p.xe(end, :);
    s.p.Ai(1, :)  = s0.p.Ai(end, :);
    s.p.Ae(1, :)  = s0.p.Ae(end, :);

    s.q.Sg(1, :)  = s0.q.Sg(end, :); 
    s.q.Sn(1, :)  = s0.q.Sn(end, :); 
    s.q.Sa(1, :)  = s0.q.Sa(end, :); 
    s.q.Hi(1, :)  = s0.q.Hi(end, :);
    s.q.He(1, :)  = s0.q.He(end, :);
    s.q.xi(1, :)  = s0.q.xi(end, :);
    s.q.xe(1, :)  = s0.q.xe(end, :);
    s.q.Ai(1, :)  = s0.q.Ai(end, :);
    s.q.Ae(1, :)  = s0.q.Ae(end, :);

end



function s = updateStateVars(s, r, dt, pars)

	tt = s.t + 1;

    [xPi, xPe, xQi, xQe] = inputCurrent(s, dt, pars);
    s.p.xi(tt, :)  = xPi;
    s.p.xe(tt, :)  = xPe;
    s.q.xi(tt, :)  = xQi;
    s.q.xe(tt, :)  = xQe;

    [dSn, dSp, dSq] = dSdt(s, r, pars);
    s.n.SPn(tt, :) = s.n.SPn(s.t, :) + dSn.Pn * dt;
    s.n.SPa(tt, :) = s.n.SPa(s.t, :) + dSn.Pa * dt;
    s.p.Sg(tt, :)  = s.p.Sg(s.t, :)  + dSp.g  * dt;
    s.p.Sn(tt, :)  = s.p.Sn(s.t, :)  + dSp.n  * dt;
    s.p.Sa(tt, :)  = s.p.Sa(s.t, :)  + dSp.a  * dt;
    s.q.Sg(tt, :)  = s.q.Sg(s.t, :)  + dSq.g  * dt;
    s.q.Sn(tt, :)  = s.q.Sn(s.t, :)  + dSq.n  * dt;
    s.q.Sa(tt, :)  = s.q.Sa(s.t, :)  + dSq.a  * dt;

    % We need to ensure S > 0 because noise can be negative
	s.n.SPn(tt, :) = heaviside(s.n.SPn(tt, :)) .* s.n.SPn(tt, :);
	s.n.SPa(tt, :) = heaviside(s.n.SPa(tt, :)) .* s.n.SPa(tt, :);
	s.p.Sg(tt, :)  = heaviside(s.p.Sg(tt, :))  .* s.p.Sg(tt, :);
	s.p.Sn(tt, :)  = heaviside(s.p.Sn(tt, :))  .* s.p.Sn(tt, :);
	s.p.Sa(tt, :)  = heaviside(s.p.Sa(tt, :))  .* s.p.Sa(tt, :);
	s.q.Sg(tt, :)  = heaviside(s.q.Sg(tt, :))  .* s.q.Sg(tt, :);
	s.q.Sn(tt, :)  = heaviside(s.q.Sn(tt, :))  .* s.q.Sn(tt, :);
	s.q.Sa(tt, :)  = heaviside(s.q.Sa(tt, :))  .* s.q.Sa(tt, :);

    [dAPi, dAPe, dAQi, dAQe] = dAdt(s, pars);
    s.p.Ai(tt, :) = s.p.Ai(s.t, :) + dAPi * dt;
    s.p.Ae(tt, :) = s.p.Ae(s.t, :) + dAPe * dt;
    s.q.Ai(tt, :) = s.q.Ai(s.t, :) + dAQi * dt;
    s.q.Ae(tt, :) = s.q.Ae(s.t, :) + dAQe * dt;

    [dHPi, dHPe, dHQi, dHQe] = dHdt(s, pars);
    s.p.Hi(tt, :) = s.p.Hi(s.t, :) + dHPi * dt;
    s.p.He(tt, :) = s.p.He(s.t, :) + dHPe * dt; 
    s.q.Hi(tt, :) = s.q.Hi(s.t, :) + dHQi * dt;
    s.q.He(tt, :) = s.q.He(s.t, :) + dHQe * dt; 

    s.t = tt;

end



function [xPi, xPe, xQi, xQe] = inputCurrent(s, dt, pars)

    J = 1:pars.N;
    I = 1:pars.N;

    % Delays
    minDelay = pars.selfDelay / dt;
    maxDelay = pars.maxDelay  / dt;
    delay = round(linspace(minDelay, maxDelay, pars.N));
    tij = s.t - delay(abs(repmat(I',size(J)) - repmat(J',size(I))') + 1);
    tij = max(1, tij);
    indPS  = sub2ind(size(s.p.Sn), tij, repmat(J, [pars.N, 1]));
    indQS  = max(1, round(s.t - minDelay));
    
    % Pitch inputs
    % Thalamic-to-pitch
    thalP = pars.Jtn * s.n.SPn(s.t, :) + pars.Jta  * s.n.SPa(s.t, :);
    % Pitch-to-pitch 
    xPei = sum(pars.Jei * pars.Cei .* s.p.Sn(indPS) + ...
               pars.Jai * pars.Cei .* s.p.Sa(indPS), 2)';
    xPii = sum(pars.Jii * pars.Cii .* s.p.Sg(indPS), 2)';
    xPee = sum(pars.Jee * pars.Cee .* s.p.Sn(indPS) + ...
               pars.Jae * pars.Cee .* s.p.Sa(indPS), 2)';
    xPie = sum(pars.Jie * pars.Cie .* s.p.Sg(indPS), 2)';
    % Q-to-pitch inputs
    zPee = pars.Jqpen * s.q.Sn(indQS, :) + pars.Jqpea * s.q.Sa(indQS, :);
    zPei = pars.Jqpin * s.q.Sn(indQS, :) + pars.Jqpia * s.q.Sa(indQS, :);
    zPie = pars.Jqpeg * s.q.Sg(indQS, :);
    zPii = pars.Jqpig * s.q.Sg(indQS, :);

    % Q inputs 
    % Pitch-to-Q
    zQee = pars.Jpqen * s.p.Sn(indQS, :) + pars.Jpqea * s.p.Sa(indQS, :);
    zQei = pars.Jpqin * s.p.Sn(indQS, :) + pars.Jpqia * s.p.Sa(indQS, :);
    zQie = pars.Jpqeg * s.p.Sg(indQS, :);
    zQii = pars.Jpqig * s.p.Sg(indQS, :);
    % Q-to-Q (same AMPA conductivities as in pitch-to-pitch)
    xQei = sum(pars.Jqei * pars.Cqei .* s.q.Sn(indPS) + ...
               pars.Jqai * pars.Cqei .* s.q.Sa(indPS), 2)';
    xQii = sum(pars.Jqii * pars.Cqii .* s.q.Sg(indPS), 2)';
    xQee = sum(pars.Jqee * pars.Cqee .* s.q.Sn(indPS) + ...
               pars.Jqae * pars.Cqee .* s.q.Sa(indPS), 2)';
    xQie = sum(pars.Jqie * pars.Cqie .* s.q.Sg(indPS), 2)';

    % Population I0s
    extE = repmat(pars.I0e, size(xPee));
    extI = repmat(pars.I0i, size(xPii));

    % Net inputs
    xPi = extI + xPei - xPii + zPei - zPii - s.p.Ai(s.t, :);
    xPe = extE + xPee - xPie + zPee - zPie - s.p.Ae(s.t, :) + thalP;

    xQi = extI + xQei - xQii + zQei - zQii - s.q.Ai(s.t, :) + pars.QI0i;
    xQe = extE + xQee - xQie + zQee - zQie - s.q.Ae(s.t, :) + pars.QI0e;

end



function [dSn, dSp, dSq] = dSdt(s, r, pars)

    n = pars.sigma * randn([1, 13]);

    leakPI = -s.p.Sg(s.t, :)  / pars.tauGABA;
    leakPE = -s.p.Sn(s.t, :)  / pars.tauNMDA;
    leakPA = -s.p.Sa(s.t, :)  / pars.tauAMPA;

    dSp.g  = leakPI + (1 / 1000) * s.p.Hi(s.t, :) + n(4);
    dSp.n  = leakPE + pars.gamma * (1-s.p.Sn(s.t, :)).*s.p.He(s.t, :) + n(5);
    dSp.a  = leakPA + (1 / 1000) * s.p.He(s.t, :) + n(6);
    
    leakQI = -s.q.Sg(s.t, :) / pars.tauGABA;
    leakQE = -s.q.Sn(s.t, :) / pars.tauNMDA;
    leakQA = -s.q.Sa(s.t, :) / pars.tauAMPA;

    dSq.g  = leakQI + (1 / 1000) * s.q.Hi(s.t, :) + n(7);
    dSq.n  = leakQE + pars.gamma * (1-s.q.Sn(s.t, :)).*s.q.He(s.t, :) + n(8);
    dSq.a  = leakQA + (1 / 1000) * s.q.He(s.t, :) + n(9);

    leakNPn = -s.n.SPn(s.t, :) / pars.tauNMDA;
    leakNPa = -s.n.SPa(s.t, :) / pars.tauAMPA;

    dSn.Pn = leakNPn + pars.gamma * (1-s.n.SPn(s.t, :)).*r.A(s.t, :) + n(10);
    dSn.Pa = leakNPa + (1 / 1000) * r.A(s.t, :) + n(11);

end



function [dAPi, dAPe, dAQi, dAQe] = dAdt(s, pars)

	dAPi = -s.p.Ai(s.t, :) / pars.tauAdapt + pars.adaptStr * s.p.Hi(s.t, :);
    dAPe = -s.p.Ae(s.t, :) / pars.tauAdapt + pars.adaptStr * s.p.He(s.t, :);
    dAQi = -s.q.Ai(s.t, :) / pars.tauAdapt + pars.adaptStr * s.q.Hi(s.t, :);
    dAQe = -s.q.Ae(s.t, :) / pars.tauAdapt + pars.adaptStr * s.q.He(s.t, :);

end



function [dHPi, dHPe, dHQi, dHQe] = dHdt(s, pars)

    [phiIpe, tauPe] = phi(s.p.xe(s.t, :), s.p.He(s.t, :), 'e', pars);
    [phiIpi, tauPi] = phi(s.p.xi(s.t, :), s.p.Hi(s.t, :), 'i', pars);
    [phiIqe, tauQe] = phi(s.q.xe(s.t, :), s.q.He(s.t, :), 'e', pars);
    [phiIqi, tauQi] = phi(s.q.xi(s.t, :), s.q.Hi(s.t, :), 'i', pars);

    dHPe = (phiIpe - s.p.He(s.t, :)) ./ tauPe; 
    dHPi = (phiIpi - s.p.Hi(s.t, :)) ./ tauPi;
    dHQe = (phiIqe - s.q.He(s.t, :)) ./ tauQe; 
    dHQi = (phiIqi - s.q.Hi(s.t, :)) ./ tauQi;

end



function [phi, tau] = phi(x, H, typ, pars)

    Delta = 1; 
    a = pars.(['a', typ]);
    b = pars.(['b', typ]);
    d = pars.(['d', typ]);

    y = a * x - b;
    phi = y ./ (1 - exp(-d * y));
    
    if pars.tauEff
        dp = a * (1 ./ (y + eps) + d ./ (1 - exp(d * y))) .* phi; 
        tau = pars.tauPop * min(1, Delta * dp' ./ H')';
    else
        tau = pars.tauPop;
    end

    tau = max(1000 / pars.cortFs, tau); % tau cannot be smaller than dt

end