1
$\begingroup$

I have a large-scale biophysical neural network, and I record the amount of spikes the cells fire during a stimulus window. There's 2 conditions: one is control and in the other one the addition of GABAzine is simulated, which blocks transmission of GABA-ergic receptors. The cells I record are under the effects of inhibition by another cell type and their synapses uses GABA, so the addition of GABAzine disinhibits my recorded cells. I calculate the difference in activity between the gabazine and control conditions for a certain cell type that has 4 dendrites which can each be connected to only 1 other cell. These presynaptic cells either transmit the input stimulus or they don't; so in the end in this input stage of the network the recorded cells can belong to 1 of 5 categories: 0, 1, 2, 3 and 4 dendrites connected to the stimulus. Cells that receive more stimulus fire more when disinhibited, so there's a linear relationship between stimulus strength and disinhibitory effect.

I've noticed however that there is a trend for the cells with 3 active dendrites to be less disinhibited and want to investigate. I have 20,000 runs with the following combinatorial design:

  • 20 trials per condition (x20)
  • with and without gabazine (x2)
  • 0, 1, 2, 3, 4 active dendrites (x5)
  • Faster/slower stimulus in 20 steps (x20)
  • Coherence between inhibitory feedback and stimulus in 5 steps (x5)

In the end I take 10,000 pairs of runs with and without gabazine and count the difference in spikes fired during the stimulus window. I end up with 100 series of this difference reported across the active dendrite categories, with 100 samples per category (5 coherence levels x 20 trials). Per sample there's 8 sources of Poisson distributed background noise on each of the 4 inhibitory and excitatory synapses that converge on the modelled cell, and the stimulus is each time drawn from a custom distribution to model a physiological "air puff" stimulus on the whiskers of a mouse. All involved biophysical models like the cells and synapses are deterministic. The cells fire about 0-5 spikes in control condition and 3-9 spikes in gabazine condition. The diff usually ranges from 1-5.

Is there anything that can confirm that category 3 is not part of the linear relationship, like drawn in this amazing boxplot?

enter image description here

The red line being the linear relationship between 1, 2, 4 and the change in activity due to GABAzine across all trials, and 3 being less than the expected linear relationship across all trials.

I followed a univerisity course on statistics that can sadly be boiled down to a decision making tree of what tests to select, but I'm very interested to learn more, so aside from the name of the test that I'd need I'd also like to learn from the answer :)

$\endgroup$
8
  • $\begingroup$ Maybe multiple regression analysis? ... but to say anything more concrete and helpful we need some specifics of your case! $\endgroup$ Apr 14, 2021 at 13:52
  • $\begingroup$ What details do you need? $\endgroup$ Apr 14, 2021 at 14:10
  • $\begingroup$ I rewrote the question for clarity $\endgroup$ Apr 14, 2021 at 14:24
  • 2
    $\begingroup$ Ok, I gave that a go :) $\endgroup$ Apr 15, 2021 at 10:08
  • 1
    $\begingroup$ done :) coherence: taken as samples, so I corrected it to 100 samples per category. y-axis: the diff in spikes between control and gabazine. All trials are independent in the same way: Each trial is simulated and only the parameters and random noise mentioned in the text differ. The cells fire about 0-5 in ctrl and 3-9 in gabazine $\endgroup$ Apr 16, 2021 at 10:30

2 Answers 2

1
$\begingroup$

Simple answer:

If we just look at the plot, what you are asking for is a simple test of linearity: are the data adequately explained by a change in response (in your plot, response is the difference in counts between +/- GABAzine) that increases linearly with the number of dendrites, or do you need a non-linear relationship?

In your case, this can be done simply by setting up the activeDendriteNumber as an ordered factor: a categorical predictor that is implicitly assumed to have equal spacing between the levels. In R, that can be accomplished by specifying ordered = TRUE when you set up the factor variable; see the manual page for factor().

You then make sure that activeDendriteNumber is coded with orthogonal polynomials in the regression; that's the default in R for an ordered categorical predictor but might not be in other software. That allows for a direct evaluation of linear, quadratic, cubic (etc., up to what's allowed by the number of category levels) contributions of activeDendriteNumber to the outcome. In R, the command

simpleModel <- lm(response ~ factor(activeDendriteNumber, ordered = TRUE)),

where lm() is the R linear regression function, would do the trick. You ask for summary(simpleModel), and you will get an evaluation of the linear, quadratic, etc contributions to the response. A significant non-linear contribution would mean that you are essentially done if the results are as you have sketched.

To do tests comparing specific situations thereafter, many find the emmeans package helpful, as it provides a coherent system for further processing of many types of models. Here, you would use simpleModel as an input to emmeans; you can also use that package with the more complicated models described below.

You might get away with this analysis, but with small numbers of counts the assumptions of standard linear regression might not hold. Furthermore, this analysis hides a lot of information that might be important for figuring out what's going on: why is there something special about an activeDendriteNumber of 3? At some point you will need to answer that question. To help figure that out, you might want to get deeper into modeling the count data from individual trials.

Modeling individual trial counts:

Try modeling the actual number of counts in each trial, not just differences +/- GABAzine, as a function of the relevant predictors: not just activeDendriteNumber as modeled above as an ordered categorical predictor, but also GABAzine presence (binary predictor, FALSE/TRUE), stimulusSpeed (continuous predictor coding is probably OK with 20 levels), and coherence (seems to be another ordered categorical predictor). Data would be organized in a long format, with each row of data representing a separate trial containing the number of counts and the predictor values in place during the trial.

As these are count data with small numbers of counts, you would use a generalized linear model. A generalized linear model estimates regression coefficients for each of a set of predictors like those in standard linear regression, but then links the value of the combined linear predictor to the outcome variable in a way commensurate with the nature of the outcome variable. Unlike standard linear regression there is no closed form for estimating the regression coefficients, so that is done with maximum-likelihood estimation.

In your case, small-valued count data are inherently non-negative and should have a variance related to the number of counts. For example, with Poisson-distributed counts, the mean equals the variance and data are modeled by the Poisson regression flavor af a generalized linear model.

In R, the glm() function with a poisson family setting and its default log link between the linear predictor and outcome counts would be a place to start, in its simplest form:

countModel <- glm(counts ~ gabazine + stimulusSpeed + activeDendriteNumber + coherence, family = poisson(link = "log"))

With the log link the regression coefficients represent differences in expected log counts, thus proportional changes in counts. This simple model might be all you need.

For example, in your plot of the +/- GABAzine count differences versus activeDendriteNumber, might counts with both GABAzine = FALSSE and GABAzine = TRUE be increasing with activeDendriteNumber? If so, then the increasing disinhibition with increasing active dendrites, based on count differences in your plot, might just be the same proportional increase with GABAzine in all cases (except for the case of 3 active dendrites). Then you can focus on whether the situation with 3 active dendrites is specifically due to changes in counts with or without GABAzine.

If the above simple model is insufficient based on your understanding of your biophysical neural network, you could evaluate interaction terms between predictors in the model. For example, if the proportionality of counts depending on GABAzine is expected to change as a function of stimulusStrength you could write:

interactionCountModel <- glm(counts ~ gabazine * stimulusSpeed + activeDendriteNumber + coherence, family = poisson(link = "log"))

where "*" represents terms for each of the 2 predictors individually plus their product (interaction term). This can be done to evaluate any effects that can't be explained by simple proportional changes in counts. In your case you have the advantage that you (presumably) can go back pretty easily to your biophysical neural network to run more experiments to test predictions that come from the statistical modeling.

$\endgroup$
6
  • $\begingroup$ For the simple answer, is there a way with a test of linearity to find out which one of the categories makes it that I need a non-linear relationship? Is it valid to just test with and without category 3 and see if suddenly the data is linear and from that infer that cat 3 is the odd one out? $\endgroup$ Apr 16, 2021 at 19:15
  • $\begingroup$ Do you happen to know of any equivalent packages in Python? Probably I can find what I need in scipy and sklearn $\endgroup$ Apr 16, 2021 at 19:16
  • $\begingroup$ @RobinDeSchepper I don't use Python so I can't help with that. My (poorly informed) sense is that Python packages can have hidden "gotchas" that are less prevalent in R. A model that just removes category 3 is different from a model with all 4 categories, with a lower number of cases and thus a lower power to detect a true remaining non-linearity. Comparisons among the 4 categories after the full model are the type of thing that emmeans in R handles well. $\endgroup$
    – EdM
    Apr 16, 2021 at 19:34
  • $\begingroup$ @RobinDeSchepper I suspect that full modeling of individual-trial count data, although it might seem more complicated at first, will end up being both simpler to explain and more informative. The "simple answer" can't tell you why category 3 is special. Also, it's possible that the linear increase (except for category 3) in +/- GABAzine count differences with active dendrite number wlll end up being a relatively constant (except for category 3) proportional change in counts in the full count model, which can then be used to explore in more detail what's happening with category 3. $\endgroup$
    – EdM
    Apr 16, 2021 at 19:40
  • 1
    $\begingroup$ @RobinDeSchepper pay particular attention to the ordered-predictor setup and orthogonal-polynomial coding in Python, whether you use GLM or not. The defaults might not be the same as in R. That said, if I'm right about constant proportional changes in counts with categories (except for category 3), you might not need the ordered coding; an unordered coding with a category * GABAzine interaction term might expose the differences associated with category 3 even more explicitly. $\endgroup$
    – EdM
    Apr 16, 2021 at 19:57
1
$\begingroup$

Not sure if understood your question, but you can test if the relationship is linear by comparing two models using for example an F test

  1. model has only the assumed linear effect y ~ intercept + beta*numberof dendrites
  2. model has indicators for groups of dendrites y ~ intercept + beta1*1dentrite + beta2*2dendrite + beta3*3dendrite

or as an alternative you can check if it is only the third group that does not follow the linear pattern by comparing model 1 or 2 to

  1. model y ~ intercept + beta1*numberofdendrites + beta2*3dendrite
$\endgroup$
1
  • $\begingroup$ Thanks, I'll try that! How do I determine the betas used in your examples? I don't think I fully understand the F-test inputs and outputs ;p $\endgroup$ Apr 16, 2021 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.