A critical roadblock to lignocellulosic biofuel is the efficient degradation of crystalline fibers of cellulose to glucose. It is caused by the unusually high thermal and mechanical stability of cellulose. The redundancy in hydrogen bonding (H-bonding) pattern and the intertwinement of intra- and inter-molecular H-bonds ensure this high stability. We have performed computations both at atomistic and coarse-grained levels to investigate the thermal responses of H-bonding networks of cellulose. We will discuss our use of all-atom simulations to understand the molecular aspects that lead to cellulose adopting its different crystal phases, and coarse-grained modeling to understand the bulk properties of microfibrils composed of those phases. All atom replica exchange molecular dynamics simulations have been used to examine (i) conformational preference of different lengths of single cellulose chains and, (ii) aggregation propensities of multiple cellulose chains. We present results that reveal the flexibility and other thermodynamic and mechanical properties. Conformational biases upon assembly are captured and compared with those of soluble cellulose chains. In the coarse-grained approach, we have constructed a statistical mechanical model at the resolution of explicit H-bonds that takes into account both intra-chain and inter-chain H-bonds in naturally occurring cellulose crystals. This model captures the plasticity of the H-bonding network in cellulose due to frustration and redundancy in available H-bonds. Furthermore, instead of only one stable H-bond pattern, different H-bonding patterns dominate at different temperatures till the disassembly at very high temperature. It provides useful clues on rational procedure for the efficient degradation.