I know that this is late but I stumbled upon this question myself too and I didn't find any answers. Fortunately, I thought about it for a while and I came up with a satisfying conclusion.
By examining the dynamic routing algorithm, one can clearly see that the logits $b_{ij}$ must always be set to zero (whether training the model or not), and then they are updated on each iteration according to the prediction $\hat{u}_{j|i}$, and in turn, update $c_{ij}$. What this means is that all $c_{ij}$ are not parameters for the model but rather values determined by the dynamic routing algorithm, hence the name "dynamic". In other words, $c_{ij}$ are not trainable parameters. Similar to weighted sums and activations of ordinary deep neural networks, you need to store these values during the feed forward pass to substitute them later on in the calculations of the backpropagation algorithm.
That said, the transformation matrices $W_{ij}$ are the trainable parameters and what you need to do is calculate their derivatives. To do that you need to substitute the equations $\hat{u}_{j|i} = W_{ij} u_{i}$ and $s_j = \sum\limits_{i} c_{ij} \hat{u}_{j|i} $ in the activation function $v_j =\frac{\lVert s_j \rVert^{2}}{1+ \lVert s_j \rVert^{2}} \frac{s_j}{ \lVert s_j \rVert}$, i.e., calculate $ s_j = \sum\limits_{i} c_{ij} W_{ij} u_{i} $ and use that to calculate $\lVert s_j \rVert$ using the euclidean norm, then plug these in $v_j$. The final formula will be daunting but at least it will be in terms of the weights, the values $c_{ij}$ which should be treated as constants, and $u_i$.
To make things a bit easier you can use this trick: set $z_{ij} = W_{ij} u_{i}$ and then the derivations should be nicer in terms of $z_{ij}$.