If I correctly understand the goal here, it sounds as though you may need to restructure your model and data slightly.
First the model:
EEG ~ 1 + age*timepoint + (1|id) + (1 + age*timepoint|node)
where node is a label that indicates which sensor position the EEG value maps onto.
Now the data you might feed this model then needs to be in "stacked" or sometimes referred to as "long" format. Each row would represent the EEG value for individual $i$ derived from the $j$th sensor.
In this scenario, I imagine you might be primarily interested in the credibility interval around the random slopes specified by the (... + age*timepoint|node) portion of the model.
For a more structural equation-y way to use brms and achieve (again what I understand) the goal to be, you could create a series of equations as follows:
form_EEG_1 <- bf(
EEG_1 ~ 1 + age*timepoint + (1|id)
) + gaussian() #assuming this is the appropriate family for EEG values
form_EEG_2 <- bf(
EEG_2 ~ 1 + age*timepoint + (1|id)
) + gaussian()
.
.
.
form_EEG_6 <- bf(
EEG_6 ~ 1 + age*timepoint + (1|id)
) + gaussian()
and then in the brm function you would add the formulae together as follows:
fit_EEG <- brm(formula = form_EEG_1 + form_EEG_2 + ... form_EEG_6 + set_rescor(TRUE),
[other brm control arguments]
)
The downside with this second approach to your model is that you would not be able to include residual correlations among your DVs if you cannot assume a multivariate normal distribution for your EEG values. So depending on what is a reasonable prior for these scores, this may or may not be an option for you.
Note these approaches are not equivalent, and they handle the dependency in your data in slightly different ways so be sure to select an approach that best maps on to your analytic goals and your data properties.